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Summary 

A  study  of  the  full-potential  modeling  of  a  blade-vortex  interaction  was  made.  A  primary  goal  of 
this  study  was  to  investigate  the  effectiveness  of  the  various  methods  of  modeling  the  vortex.  The 
problem  was  studied  within  the  context  of  a  two-dimensional  model  problem,  which  represents  one  of 
the  limiting  types  of  blade- vortex  interactions  (BVIs).  The  model  problem  restricts  the  interaction 
to  that  of  an  infinite  wing  with  an  infinite  line  vortex  moving  parallel  to  its  leading  edge.  This 
problem  provides  a  convenient  testing  ground  for  the  various  methods  of  modeling  the  vortex  while 
retaining  the  essential  physics  of  the  full  three-dimensional  interaction.  The  flow  field  is  assumed  to 
be  inviscid,  irrotational,  unsteady,  and,  in  general,  transonic. 

A  full- potential  algorithm  specifically  tailored  to  solve  BVI  was  developed  to  solve  this  problem. 
The  algorithm  makes  use  of  the  unsteady  mass  conservation  and  Bernoulli  equations  to  form  a  full- 
potential  model  of  the  flow  field.  The  system  of  equations  is  reduced  to  one  equation  by  using  a 
Taylor-series  expansion  of  the  temporal  derivative  of  the  density  term  in  the  conservation  equation. 
The  spatial  derivatives  are  recast  in  “delta”  form,  with  the  density  written  at  the  previous  time  step. 
The  stability  of  the  algorithm  in  transonic  flow  is  assured  through  the  use  of  upwind  biasing  of  the 
density  in  the  flux  terms.  The  flux  metrics  are  computed  by  the  consistent  metric  method,  which 
has  been  found  to  be  superior  to  the  so-called  free-stream  subtraction  method  that  has  difficulties 
with  grid  singularities.  The  equation  is  approximately  factored  into  convenient  geometric  parts  in 
order  to  reduce  the  matrix  to  a  compact  form.  A  tridiagonal  matrix  inversion  is  used  to  solve  for 
the  updated  potential  solution.  The  model  has  the  capability  to  predict  the  steady  and  unsteady 
flow  about  an  airfoil  under  subcritical  and  transonic  flow  conditions.  Comparisons  of  the  results 
predicted  are  made  with  those  presented  by  other  researchers  and  with  experimental  data.  The 
comparisons  indicate  that  the  algorithm  is  able  to  predict  basic  unsteady  transonic  flow  about  an 
airfoil. 

The  basic  algorithm  has  been  modified  to  include  the  effect  of  a  vortex  passing  near  the  airfoil. 
Four  different  methods  of  modeling  of  the  vortex  were  used: 

1.  The  angle- of- at  tack  method 

2.  The  lifting-surface  method 

3.  The  branch-cut  method 

4.  The  split- potential  method 

The  angle-of-attack  method  uses  the  velocity  field  of  a  point  vortex  to  compute  a  vortex-induced 
velocity  at  the  airfoil  quarter-chord.  This  velocity  is  then  used  to  compute  an  effective  angle 
of  attack  of  the  airfoil.  This  method  is  identical  to  techniques  which  are  currently  in  use  in 
comprehensive  helicopter  rotor  analyses.  The  lifting-surface  method  is  an  extension  of  the  angle-of- 
attack  method  in  which  the  vortex -induced  velocity  is  a  function  of  chordwise  distance  on  the  airfoil 
surface.  The  branch-cut  method  is  a  flow-field  vortex  representation  that  makes  use  of  a  surface  of 
potential  discontinuity,  the  edge  of  which  constitutes  the  vortex  location.  The  effect  of  the  vortex  is 
implemented  by  imposing  special  differencing  methods  along  the  cut.  In  the  split- potential  method, 
the  velocity  field  is  split  between  a  known  field  (induced  by  the  vortex)  and  an  unknown  perturbation 
field  caused  by  the  airfoil. 

A  side-by-side  comparison  of  the  four  models  was  conducted.  These  comparisons  include 
comparing  generated  velocity  fields,  a  subcritical  interaction,  and  a  critical  interaction.  The 
subcritical  and  critical  interactions  are  compared  with  experimentally  generated  results. 

The  split-potential  model  was  used  to  make  a  survey  of  some  of  the  more  critical  parameters 
which  affect  the  BVI.  The  survey  studies  general  flow  parameters  such  as  free-stream  Mach  number 
and  airfoil  angle  of  attack,  and  vortex  parameters  such  as  strength,  core  size,  and  miss  distance.  The 
results  were  computed  at  subcritical  and  supercritical  free-stream  Mach  numbers.  For  the  vortex 
parameters,  the  free-stream  Mach  number  was  chosen  to  be  just  subcritical  in  order  to  study  the 
effect  of  the  vortex  on  the  formation  of  critical  flow  on  the  airfoil. 


XI 


1.  Introduction 


1.1.  Physical  Problem 

Helicopter  rotors  operating  in  high-speed  flight 
encounter  a  number  of  important  aerodynamic  phe¬ 
nomena.  Two  major  features  that  dominate  the  flow 
on  the  advancing  side  of  the  rotor  disc  exist.  The 
first  key  feature  is  the  presence  of  transonic  flow  con¬ 
ditions.  Transonic  flow  imposes  major  limitations  on 
the  high-speed  performance  of  the  rotor.  These  lim¬ 
itations  manifest  themselves  in  high  vibration  levels, 
power  divergence,  noise,  and  component  fatigue.  The 
second  key  feature  is  the  presence  of  a  large  vortex 
system  near  the  rotor.  The  vortex  system  is  com¬ 
posed  of  a  series  of  helical  vortex  filaments  generated 
at  the  tip  of  the  blade.  The  following  blade,  which 
may  be  experiencing  transonic  flow,  frequently  inter¬ 
acts  with  these  vortices. 

The  interaction  between  a  rotor  blade  and  the 
vortices  from  a  preceding  blade  can  have  a  large  im¬ 
pact  on  the  blade  aerodynamic  environment.  These 
blade-vortex  interactions  (BVIs)  cause  large  changes 
in  local  pressure  which  can  occur  over  short  periods. 
The  pressure  changes  exist  within  a  flow  field  which 
is,  in  general,  transonic,  unsteady,  viscous,  and  three- 
dimensional.  The  vortex  passage,  therefore,  acts  to 
modify  an  already  complicated  flow  field. 

A  rotor  interacts  with  a  vortex  under  a  wide 
range  of  relative  orientations.  However,  the  essential 
physics  can  be  illustrated  by  considering  a  rectangu¬ 
lar  blade  of  infinite  aspect  ratio  interacting  with  an 
infinite  line  vortex  at  an  angle  0 .  Johnson  (ref.  1)  has 
shown  that  this  problem  is  steady  in  a  coordinate  sys¬ 
tem  whose  origin  travels  with  the  intersection  of  the 
blade  centerline  and  the  projection  of  the  free  vortex 
on  the  blade.  (See  fig.  1.)  The  steady  coordinate 
system  is 

x'  =  X  I 

y  =  y  -  Moo  cot  6^  J  ^  ^ 

The  speed  at  which  the  origin  travels  is  a  function 
of  the  angle  9.  When  9  =  7r/2,  the  vortex  is  per¬ 
pendicular  to  the  blade  and  the  speed  of  the  inter¬ 
action  point  is  zero.  (See  fig.  2.)  For  increasing 
values  of  9,  the  speed  of  the  interaction  point  in¬ 
creases,  but  the  problem  remains  steady.  For  9  =  tt, 
there  is  no  spanwise  flow  dependence,  and  the  prob¬ 
lem  is  now  two-dimensional.  (See  fig.  3).  How¬ 
ever,  the  cost  of  this  two-dimensional  simplification 
is  that  this  problem  is  now  intrinsically  unsteady  be¬ 
cause  the  speed  of  the  interaction  point  is  infinite. 
The  blade-vortex  interaction  may  then  be  classified 
by  the  two  limiting  conditions  defined  by  d  =  7r/2 
and  9  =  IT.  The  first  condition  (d  =  7r/2)  may  be 


Figure  1.  Interaction  of  infinite  aspect  ratio  blade  with  infi¬ 
nite  line  vortex. 


r 


Figure  2.  Low-speed  interaction  between  rotor  and  vortex. 


0  =  tt;  =  oo 


Figure  3.  High-speed  interaction  between  rotor  and  vortex. 

called  a  low-speed  interaction  (LSI)  because  this  is  a 
steady  problem  even  in  the  original  coordinate  sys¬ 
tem.  The  second  condition  (9  =  tt)  may  be  called 
a  high-speed  interaction  (HSI)  because  this  is  an 
unsteady  problem  even  in  the  transformed  coordinate 


system.  Both  LSI  andHSI  represent  real  interactions 
which  can  have  significant  effects  on  the  rotor  aero¬ 
dynamics.  LSI,  for  instance,  is  the  principal  type  of 
interaction  which  occurs  during  hovering  flight.  LSI 
affects  rotor  power  and  low  harmonic  loading.  HSI 
occurs  during  high-speed  flight  and  descents  and  af¬ 
fects  noise,  vibrations,  and  the  higher  harmonics  of 
loading.  Furthermore,  HSI  contains  all  the  physics 
of  LSI;  therefore,  the  capability  to  solve  for  HSI  con¬ 
tains  the  ability  to  solve  LSI.  The  solution  of  HSI  is 
the  main  study  of  this  report. 

1.2.  Physical  Model 

The  solution  of  HSI  requires  the  computation  of 
the  time-varying  surface  pressures  during  the  vor¬ 
tex  passage.  Because  the  angle  between  the  vor¬ 
tex  axis  and  the  blade  is  zero,  no  spanwise  flows 
are  induced  and  the  flow  can  be  assumed  to  be  two- 
dimensional.  The  distance  between  the  airfoil  and 
the  vortex  is  assumed  to  be  large  enough  to  assure 
a  basically  in  viscid  flow;  that  is,  the  vortex  does  not 
distort  the  airfoil  boundary  layer.  Any  shock  waves 
present  are  assumed  to  be  weak  and  not  a  source  of 
rotational  flow.  With  these  assumptions,  the  aerody¬ 
namic  problem  can  be  modeled  by  assuming  a  two- 
dimensional  potential  flow  field.  The  mass  conserva¬ 
tion  equation  for  such  a  flow  field  is 

I +  »  (12) 

Since  the  flow  is  isentropic,  the  fluid  density  can  be 
determined  from  the  Bernoulli  equation 


P  = 


(1.3) 


These  two  equations,  solved  together,  constitute  a 
full- potential  model  of  the  flow  field.  A  solution  for 
these  equations  can  be  accomplished  by  using  a  finite- 
difference  algorithm.  Such  an  algorithm,  originally 
developed  by  Steger  and  Caradonna  (ref.  2)  is  pre¬ 
sented  in  section  2.  This  scheme  is  then  modified  to 
include  various  vortex  models. 


1.3.  Vortex  Models 


A  primary  aim  of  this  report  is  to  study  the  effect 
of  various  vortex  models  on  the  solution  of  HSI  within 
the  framework  of  a  2-D  potential  finite-difference 
algorithm.  All  these  models  are  candidates  for  use 


in  3-D  methods.  Four  different  methods  have  been 
used  to  model  the  vortex: 

1.  The  first  model  approximates  the  vortex  effect 
as  a  change  in  airfoil  angle  of  attack.  The 
velocity  field  of  a  point  vortex  law  is  used 
to  compute  an  induced  velocity  at  the  airfoil 
quarter  chord.  This  induced  velocity  is  then 
applied  over  the  entire  airfoil.  This  is  referred 
to  as  the  “angle-of-attack  method.” 

2.  The  second  model  is  related  to  the  first  but, 
instead  of  imposing  only  a  constant  velocity  on 
the  airfoil,  a  distributed  velocity  field  from  the 
vortex  is  imposed  on  the  airfoil  surface.  This 
is  analogous  to  a  lifting-surface  method. 

3.  The  third  model  is  to  specify  a  branch-cut 
discontinuity  in  the  potential  field.  The  vortex 
is  modeled  as  a  jump  in  potential  across  the 
branch  cut,  the  edge  of  which  represents  the 
center  of  the  vortex.  This  is  referred  to  as  the 
“branch-cut  method.” 

4.  The  fourth  method  models  the  vortex  by  ex¬ 
pressing  the  potential  as  the  sum  of  a  known 
potential  due  to  the  vortex  and  an  unknown 
potential  due  the  airfoil.  This  is  referred  to  as 
the  “split-potential  method.” 

The  first  two  vortex  models  are  typical  of  the  linear 
integral  flow  methods  which  are  used  in  all  the  cur¬ 
rently  available  rotor- analysis  methods.  Methods  3 
and  4  are  flow  models  of  the  vortex  and  can  only  be 
used  in  finite- difference  methods.  An  important  as¬ 
pect  of  this  work  is  to  determine  whether  this  more 
elaborate  modeling  is  necessary. 

1.4.  Historical  Background 

The  problem  of  blade-vortex  interactions  is  cen¬ 
tral  to  helicopter  aerodynamics  because  the  interac¬ 
tion  of  the  rotor  and  its  vortex  system  can  have  a 
large  effect  on  the  aerodynamic  environment  of  the 
blade.  The  four  vortex  models  discussed  in  section 
1.3  reflect  the  level  of  sophistication  of  the  global 
theories  within  which  they  were  developed.  To  ap¬ 
preciate  fully  the  various  vortex  models,  reviewing 
the  basics  of  the  global  rotor  computations  is  useful. 

Because  of  the  geometric  complexity  of  the  vor¬ 
tex  system,  early  analysts  (e.g.,  Glauert  in  1926 
(ref.  3))  treated  the  wake  influence  on  the  rotor  by 
using  momentum  theory  and  blade  element  approxi¬ 
mations.  The  resulting  models  led  to  simple  algebraic 
equations  for  the  induced  velocities.  Computing  an 
effective  angle  of  attack  on  each  blade  segment  with 
the  induced  velocity  at  the  disc  is  possible  and  is  the 
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essence  of  the  angle- of- attack  method.  The  result¬ 
ing  angle  of  attack  is  then  used  in  conjunction  with 
tabulated  lift,  drag,  and  moment  data  to  compute 
rotor  performance.  This  approach  has  proven  to  be 
useful  in  computing  total  aircraft  performance  but 
is  unable  to  predict  accurately  the  details  of  the  ro¬ 
tor  aerodynamic  environment.  Accurate  prediction 
of  these  details  could  not  be  performed  until  the  ad¬ 
vent  of  high-speed  computers  (circa  1960). 

Piziali  and  DuWalt  (ref.  4)  proposed  the  first 
practical  method  for  studying  the  details  of  the  rotor- 
wake  interaction.  The  complete  rotor  wake  was 
modeled  with  a  series  of  such  line  segments.  The 
velocities  induced  by  each  segment  at  the  rotor  were 
combined  to  produce  an  effective  angle  of  attack.  The 
solution  then  proceeds  as  before.  This  improved  the 
earlier  model  by  allowing  for  individual  blade- vortex 
encounters  to  be  studied.  Isay  (ref.  5)  presented  a 
more  general  solution  of  the  induced  velocity  for  a 
spiral  wake  model. 

Kocurek  (ref.  6)  presented  an  extension  of  this 
method  (for  hover  only)  in  which  the  blade  was 
treated  as  a  lifting  surface.  The  induced  velocity  was 
computed  over  the  entire  surface  of  the  blade  not 
just  at  one  point.  These  velocities  were  then  used  to 
compute  local  pressure  and  lift.  Tabulated  data  were 
used  to  provide  the  drag  and  pitching  moment  with 
the  solution  proceeding  as  before. 

These  approaches  can  be  broadly  classified  as 
blade-element  integral  methods.  Currently  they  are 
the  most  popular  methods  of  rotor  analysis  in  use 
(especially  the  angle-of-attack  method).  Numerous 
investigators  have  improved  on  the  basic  model  cul¬ 
minating  in  the  effort  by  Johnson  (ref.  7).  Despite  an 
impressive  versatility,  however,  the  blade-element  in¬ 
tegral  methods  have  some  serious  shortcomings.  Pri¬ 
marily  the  deficiencies  of  these  methods  he  in  the 
use  of  tabulated  airfoil  data  to  provide  aerodynamic 
forces  and  in  the  related  assumption  of  local  two- 
dimensional  flow.  Furthermore,  unsteady  aerody¬ 
namic  effects  are  modeled  with  quasi-steady  approx¬ 
imations  which  are  incapable  of  modeling  the  truly 
unsteady  phenomena  of  transonic  flow. 

During  the  time  period  1960-1980,  the  field 
of  computational  fluid  dynamics  underwent  rapid 
growth.  Computer  speeds  increased  to  the  point 
where  it  became  possible  to  use  finite-difference 
methods  to  compute  simple  rotor  flows.  Various  in¬ 
vestigators  (e.g.,  refs.  8  and  9)  addressed  the  lim¬ 
itations  of  the  integral  methods  by  using  finite- 
difference  methods  to  compute  rotor  aerodynamics, 
typically  with  hybrid  methods.  The  hybrid  methods 


use  finite-difference  techniques  to  solve  a  limited  part 
of  the  flow  field  and  a  linear  integral  method  to  pro¬ 
vide  the  global  solution.  This  method  provides  the 
capability  to  compute  the  entire  transonic  nonlinear 
flow  field  near  the  rotor.  An  essential  difference  be¬ 
tween  the  nonlinear  and  the  linear  integral  parts  of 
the  solution  is  that  the  linear  solutions  depend  only 
on  the  blade  surface  and  shear  layer  conditions  be¬ 
cause  the  speed  of  sound  is  assumed  to  be  constant. 
In  contrast  the  nonlinear  solution  depends  on  flow 
conditions  in  the  entire  flow  field.  For  rotors,  this 
field  dependence  is  especially  important  because  the 
field  is  frequently  occupied  with  vortices  from  pre¬ 
vious  blades.  Therefore,  an  important  part  of  the 
development  of  rotor  finite-difference  schemes  is  the 
means  of  specifying  vortices. 

Caradonna,  Tung,  and  Desopper  (ref.  8)  devel¬ 
oped  the  first  finite-difference  scheme  that  included 
vortices  in  the  flow  field.  They  solved  a  high-tip- 
speed  hover  problem  in  which  the  vortices  were  spec¬ 
ified  as  edges  of  potential  discontinuities  (branch-cut 
method).  This  scheme  produced  good  comparisons 
with  pressure  data.  Interestingly,  reference  8  also  re¬ 
ported  an  inability  to  obtain  a  good  solution  when 
the  effect  of  the  vortex  was  included  only  by  a  blade 
surface  inflow  specification  (angle-of-attack  method). 
Strawn  and  Caradonna  (ref.  9)  solved  a  similar  prob¬ 
lem  by  using  a  split-potential  model.  Their  method  is 
a  modified  version  of  a  full- potential  algorithm  devel¬ 
oped  by  Bridgeman,  Steger,  and  Caradonna  (ref.  10). 

To  date,  the  forward  flight  computations  have  re¬ 
lied  on  vortex-induced  surface  inflow  boundary  con¬ 
ditions  (angle-of-attack  method)  and  have  been  fairly 
successful  at  high  advance  ratios  where  the  induced 
flow  is  a  small  percentage  of  the  total  inflow.  Never¬ 
theless,  there  remains  a  serious  question  of  how  best 
to  introduce  moving  vortices  into  a  computational 
grid  and  thereby  predict  their  effect.  The  solution  of 
the  2-D  HSI  is  a  convenient  testing  ground  for  the 
vortex  modeling  schemes  that  are  required  for  the 
full  3-D  problem. 

A  number  of  investigators  (refs.  11  to  17)  have 
studied  the  2-D  HSI  problem  by  using  finite- 
difference  methods.  These  investigators  have  been 
primarily  interested  in  acoustic  effects  and  have  used 
modified  versions  of  earlier  algorithms  and  vortex 
models.  The  problem  of  vortex  specification  has  not 
been  a  primary  aim  of  these  studies.  George  and 
Chang  (ref.  11)  modeled  the  vortex  with  the  angle- 
of-attack  method  to  investigate  the  effects  of  blade- 
vortex  interactions  on  noise.  Later  they  extended 
their  methodology  and  results  to  reflect  the  results 
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from  the  branch-cut  method  (ref.  12).  McCroskey 
and  Goorjian  (ref.  13)  introduced  the  split-potential 
method,  which  was  first  proposed  by  Steinhoff.  (See 
ref.  14.)  Computations  with  Euler  and  thin-layer 
Navier-Stokes  algorithms  have  been  presented  in  ref¬ 
erences  15  and  16.  Sankar  and  Malone  (ref.  17)  pre¬ 
sented  a  full-potential  solution  by  using  a  strongly 
implicit  procedure.  All  these  methods  produce  re¬ 
sults  which  are  characteristically  similar  and,  to 
a  limited  extent,  show  good  correlation  with  each 
other,  especially  those  methods  that  employ  the  split- 
potential  vortex  model. 


full-potential  flow  model;  this  will  be  accomplished 
for  the  first  time  in  the  present  work.  Details  of  the 
vortex  modeling  are  presented  in  section  3. 

The  third  tasks  is  a  unique  side-by-side  com¬ 
parison  between  the  models.  Comparisons  with 
experimental  data  for  sub  critical  and  critical  flow 
conditions  are  also  made.  These  comparisons  high¬ 
light  the  merits  of  the  various  models.  This  work 
provides  the  basis  for  a  more  systematic  approach 
to  three-dimensional  computations  of  blade-vortex 
interaction.  The  results  of  these  comparisons  are  pre¬ 
sented  in  section  4. 


1.5.  Purpose  of  Current  Research 

Of  all  these  methods,  the  full-potential  approach 
is  probably  the  best  suited  to  rotor  computations 
because  it  is  geometrically  general  and  is  much  faster 
than  the  Euler  and  Navier-Stokes  methods.  The 
purpose  of  the  current  research  is  to  explore  fully 
the  BVI  phenomena  within  the  context  of  the  full- 
potential  algorithm.  This  exploration  requires  four 
major  tasks. 

The  first  task  is  to  develop  a  full-potential  algo¬ 
rithm  which  is  specifically  tailored  to  solve  the  BVI 
problem.  Improvement  to  previous  methods  include 
(1)  special  boundary  conditions  to  increase  flexibil¬ 
ity  in  modeling  airfoils,  (2)  a  “full-potential”  mesh, 
(3)  allowance  for  a  variable  time  step  in  the  unsteady 
solution,  and  (4)  an  improved  method  of  comput¬ 
ing  the  flux  metric.  Details  of  these  improvements 
are  presented  in  section  2,  which  includes  a  complete 
derivation  of  the  full-potential  algorithm. 

The  second  task  is  to  implement  the  various  meth¬ 
ods  of  modeling  a  vortex.  This  task  can  best  be 
accomplished  within  the  framework  of  a  single  algo¬ 
rithm.  Using  a  single  algorithm  eliminates  any  ques¬ 
tion  of  differences  raised  by  issues  such  as  grid  or 
boundary  conditions.  The  present  method  provides 
a  unique  opportunity  to  accomplish  this.  The  full- 
potential  algorithm  is  modified  to  include  the  effects 
of  four  different  ways  to  compute  the  influence  of  the 
vortex.  These  modifications  include 

1.  Modifications  to  airfoil  boundary  conditions 
for  the  angle- of-at  tack  and  lifting-surface 
methods 

2.  Special  internal  boundary  conditions  to  imple¬ 
ment  the  branch-cut  method 


The  fourth  task  is  to  make  a  parametric  study  of 
the  BVI  problem,  which  will  be  the  first  complete 
study  presented  in  the  literature.  The  parameters 
studied  fall  into  two  categories:  (1)  flow-field  param¬ 
eters  and  (2)  vortex  parameters.  The  results  of  these 
parametric  studies  are  presented  in  section  4. 


2.  Full-Potential  Algorithm 


The  aerodynamic  problem  of  HSI  will  be  modeled 
with  a  potential  flow-field  assumption.  Under  this 
assumption  the  basic  equations  of  fluid  dynamics 
(mass,  momentum,  energy,  and  equation  of  state) 
are  reduced  to  a  system  of  two  equations  with  two 
unknowns:  the  mass  conservation  equation, 

dp  d  d  ,  .  ,  , 

and  the  Bernoulli  equation. 


P  = 


1  +  - — -  ( -  $2 


(2.2) 


In  these  equations,  all  velocities  are  normalized  by 
Goo]  distances,  by  the  airfoil  chord;  time,  by  c/aoo] 
and  density,  by  its  free-stream  value. 


2.1.  Conservative  Formulation  of 
Transformed  Equation 


The  system  of  equations  (eqs.  (2.1)  and  (2.2)) 
is  transformed  to  a  computational  plane  under  the 
general  transformation 


3.  The  inclusion  of  a  split-potential  model 

The  split-potential  method  has  not  yet  been  fully 
implemented  for  unsteady  BVI  problems  with  the 


11=  ri{x,y,t)\  (2.3) 

T=t  ] 


4 


and  the  equations  must  be  conservative  in  these 
coordinates .  Equation  (2.1)  is  written  in  conservative 
(or  divergence)  form;  that  is, 

(9F  ■ 

where  F,  is  the  flux  of  the  quantity  /  {p  for  the 
potential  model).  This  generic  form  must  be  main¬ 
tained  through  the  entire  solution  process  (includ¬ 
ing  discretization)  if  mass  conservation  is  maintained. 
Viviand  (ref.  18)  has  derived  a  general  conservation 
form  for  such  a  generalized  coordinate  system  which 
transforms  equation  (2.1)  into 


dr 


(2.5) 


where  J„  is  the  Jacobian  of  the  transformation 
(eqs.  (2.3)).  Under  this  transformation,  the  Bernoulli 
equation  becomes 


P  = 


2^r 


(2.6) 


where  U  and  V  are  the  contravariant  velocity  vec¬ 
tors,  with  U  being  the  velocity  perpendicular  to  the 
r]  direction  and  V  being  the  velocity  perpendicular  to 
the  ^  direction.  In  general  these  velocities  are  defined 
as 

U  =  ] 

}  (2.7) 

V  =r/^+ A2$^+ 

where  Ap  A2,  and  A3  are  metric  terms  defined  as 


with  the  current  method.  The  streamlines  and  po¬ 
tentiallines  which  surround  an  airfoil  in  incompress¬ 
ible  flow  form  such  a  grid.  This  type  of  grid  may 
be  computed  by  means  of  a  complex  mapping  solu¬ 
tion.  The  Joukowski  airfoil  transformation  is  used  in 
the  present  method  because  it  provides  a  convenient 
closed-form  solution.  The  grid  is  generated  with  the 
following  steps: 

1.  Produce  a  satisfactory  stretched  Cartesian 
grid  using  any  method 

2.  Use  the  (^,?/)  coordinates  along  the  front  face 
of  the  grid  to  integrate  ^  to  the  aft  face  of  the 
grid  (this  solves  for  the  streamlines  around  a 
circle) 

3.  Transform  the  circle  solution  by  using  the 
Joukowski  transformation  to  produce  the  air¬ 
foil  solution 

4.  Select  an  appropriate  distribution  of  points 
along  the  airfoil  “streamline” 

5.  Interpolate  $(^,r/)  to  And  the  potential  at 
each  of  these  points 

6.  Find  the  location  of  each  of  the  “off  air¬ 
foil”  streamlines  which  have  matching  poten¬ 
tial  values 

7.  Form  the  grid  with  the  resulting  set  of  points 

The  details  of  the  development  are  presented  in 
appendix  A.  Since  the  grid  is  orthogonal,  the  metric 
term  A2  is  identically  zero.  Furthermore,  since  the 
grid  is  steady,  the  contravariant  velocities  become 


Ai  =  Ve  •  Ve  =  d  +  1 

A2=V^-Vr]  =^xr]x  +  iyr]y  \  (2.8) 

A3  =Vr]  -Vt]  =  +  r]y  J 


U  = 

V  =  A3$J 

2.3.  Metric  Derivatives 


(2.9) 


2.2.  Computational  Grid 

For  the  present  purposes,  three  characteristics 
are  useful  for  the  computational  grid  to  have.  The 
first  characteristic  is  orthogonality,  which  is  useful 
because  it  reduces  the  complexity  of  an  algorithm. 
The  second  characteristic  is  that  the  grid  lines  con¬ 
form  closely  to  the  shape  of  the  airfoil;  this  increases 
the  accuracy  of  the  solution.  The  third  character¬ 
istic  is  that  the  grid  lines  should  align  with  the  free 
stream  away  from  the  surface  in  order  to  facilitate  the 
branch-cut  method  of  modeling  the  vortex.  For  these 
reasons,  an  orthogonal  H-mesh  was  chosen  to  be  used 


The  transformation  of  the  equations  to  the  com¬ 
putational  plane  gives  rise  to  metric  terms  (i.e., 
rfy)  as  a  result  of  the  chain  rule.  In  matrix  form,  the 
chain  rule  expansion  is 


Xr  yr  1 

dx 

dr 

x^  y^  0 

dy 

= 

d^ 

yri  0  _ 

.dy. 

The  determinant  of  the  matrix  is 

Dn  =  -  Xj^y^  (2-11) 
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Cramer’s  rule  can  be  used  to  solve  for 
yjjd^  -  y^dr, 


dx  — 


Dn 

Ojj  X  0^ 


(2.12) 


«,=  O 


These  expressions  are  applied  to^  and  r],  respectively, 
to  obtain 


yy 

ivx) 

'  T  ^ 
^7] 

iVy) 

~  ^ 

J  n 

=  D- 

1 

and  X, 
( 

t  are 

(2.13) 


metrics  and  are  determined  by  the  following  simple 
central  difference  formulas  (for  convenience  = 
At]  =  1) : 


(®{)*  J 

j 

{yn)ij 


2 

(i/i+1  -  yj-l) 
2 

(®j+l  ~  l) 
2 

(r/j+1  -  r/j-i) 
2 


(2.14) 


At  the  airfoil  boundaries,  the  derivatives  in  the 
direction  are  computed  with  the  aid  of  a  pseudogrid 
line  inside  the  airfoil  contour.  This  grid  line  is 
determined  with  a  simple  linear  extrapolation  of  the 
grid  points  off  the  airfoil. 


2.4.  Boundary  Conditions 


The  four  boundary  conditions  that  are  imposed 
on  the  flow  are  associated  with  (1)  the  airfoil  surface, 
(2)  the  outer  boundary  of  the  grid,  (3)  the  aft  face  of 
the  grid,  and  (4)  the  Kutta  condition. 

Airfoil  Surface 


For  inviscid  flow,  the  surface  boundary  condition 
requires  that  the  flow  be  tangent  to  the  airfoil  sur¬ 
face.  This  requirement  can  be  met  by  setting  the 
contravariant  velocity  vector  V  to  zero.  For  a  mesh 


which  exactly  conforms  to  the  airfoil,  this  leads  to 
=  0.  One  problem  with  employing  this  boundary 
condition  is  that  a  new  mesh  must  be  generated  with 
every  new  airfoil  or  airfoil  orientation.  An  alternative 
to  computing  a  new  grid  is  to  use  a  transpiration 
rather  than  a  “no-flow”  boundary  condition.  This 
approach  uses  a  fixed  grid  which  conforms  to  some 
convenient  profile  (e.g.,  a  Joukowski  airfoil)  to  ap¬ 
proximate  the  desired  profile.  The  flow  must  there¬ 
fore  pass  through  the  grid  surface  at  an  angle  e,  which 
is  the  angular  difference  between  the  grid  surface  and 
the  actual  surface.  Figure  4  illustrates  this  relation¬ 
ship.  The  flow  normal  to  the  grid  surface  is 

$j\r  =  $s  tan  e  (2-15) 

where  is  the  flow  velocity  tangent  to  the  grid 
surface.  This  condition  is  merely  a  generalization 
of  the  usual  small-disturbance  boundary  condition. 
The  value  of  to  be  used  in  the  actual  algorithm 
remains  to  be  found.  Since  the  coordinate  system 
is  orthogonal,  the  only  difference  between  the  N  and 


(b)  Flow-through  condition  on  coordinate  which  approxi¬ 
mates  body  (generalfeation  of  small-distnrbanceboundary 
condition). 


Figure  4.  Transpiration  boundary  condition  at  airfoil  surface. 
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r]  directions  is  a  simple  one-dimensional  stretching; 
therefore, 

—  (2.16) 
r]N 

2.4-2.  Outer  Boundaries 

Along  the  outer  boundaries,  the  flow  is  required  to 
return  to  undisturbed  conditions.  Often,  computa¬ 
tional  outer  boundaries  are  so  close  to  the  airfoil  sur¬ 
face  that  this  condition  cannot  produce  an  accurate 
or  stable  result.  For  these  close-grid  boundaries,  spe¬ 
cial  nonreflective  boundary  conditions  are  imposed. 
However,  the  current  grid  boundaries  are  sufficiently 
far  away  (155  chords  horizontally  and  80  chords  ver¬ 
tically)  that  the  assumption  of  undisturbed  flow  is 
valid.  The  outer  boundary  conditions  are  set  with  a 
Dirichlet  condition 

^=M^x  (2.17) 

2.4.  3.  Aft  Faee 

Along  the  aft  face  of  the  mesh,  the  flow  is  also 
required  to  be  undisturbed.  However,  because  the 
present  method  employs  a  number  of  branch  cuts 
(lines  of  potential  discontinuity  which  model  vortic- 
ity),  the  potential  cannot  be  easily  specified  at  this 
boundary.  Instead  free-stream  conditions  are  im¬ 
posed  by  modifying  the  outgoing  flux  along  the  aft 


face  so  that  p  =  1  is  ensured.  With  the  Bernoulli 
condition,  the  following  expression  for  $ ^  is  derived 
(see  appendix  B) : 

=  —  (  Moo  -  )  (2.18) 

^  ^  Moo  J  ^  ’ 

which  is  used  in  the  flux  computation. 

2.4. 4-  Kutta  Condition 

For  the  lifting  conditions,  allowance  must  be  made 
for  a  jump  in  potential  across  a  wakelike  branch  cut 
(Kutta  condition).  This  cut  extends  from  the  air¬ 
foil  trailing  edge  to  the  aft  face  of  the  mesh  (which 
precludes  the  use  of  Dirichlet  boundary  conditions 
along  the  aft  face).  The  cut  is  aligned  with  the  mean 
chord  line  of  the  airfoil.  In  an  unsteady  flow,  the 
jump  in  potential  F  across  the  cut  must  be  convected 
downstream.  The  following  equation,  governing  the 
convection  of  vorticity  from  the  trailing  edge,  is  de¬ 
rived  by  using  the  Bernoulli  equation  and  continuity 
of  density  across  the  cut  (see  appendix  B): 

F^  +  (U)F^  =  0  (2.19) 

where  (U)  is  the  average  of  the  velocities  above  and 
below  the  branch  cut.  Equation  (2.19)  is  used  to 
determine  the  value  of  F  along  the  branch  cut. 


2.5.  Differ  ence  Equation 


Two  criteria  are  useful  in  developing  the  difference  equation.  The  first  criterion  is  that,  at  any  current  time 
step,  the  difference  equation  be  a  function  of  $  only.  This  criterion  eliminates  the  need  for  solving  a  system 
of  equations.  The  second  criterion  is  that  the  conservation  form  be  maintained.  This  criterion  ensures  an 
accurate  computation  of  shock  motion  and  strength. 

Equation  (2.5)  may  be  written  by  using  a  backward  difference  in  time  as 


n  \  /  n  \  ^ 

—  ]  -  (  —  ) 


(py_ 

VJ 


=  0 


(2.20) 


This  difference  equation  maintains  the  conservation  form,  but  the  density  still  remains  at  the  (N  +  1)  time  level; 
this  can  be  corrected  by  expanding  the  density  in  a  Taylor  series  expansion  as 


#-l-l 


-  $^) 


(2.21) 


The  operator  is  derived  from  the  Bernoulli  equation  and  is  expressed  in  the  following  conventional  form: 

dp\^  o  o  f  d  ua  va\^ 


(2.22) 
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A  complete  derivation  of  equation  (2.22)  is  provided  in  appendix  C  for  students  who  may  be  unfamiliar  with 
the  form.  The  Taylor  series  must  be  taken  on  both  the  and  terms  in  order  to  maintain  conservation 

form.  Application  of  the  appropriate  time  derivatives  dr  leads  to  the  following  dijferential  equation: 


#+l 


N 


N 


liv 


N-1 


N-1 


(2.23) 


The  terms  in  brackets  in  equation  (2.23)  are  seen  to  have  a  form  similar  to  the  time  derivatives  in  the  following 
non  conservative  full- potential  equation: 


$^^+2U$^^+2V$^^  =  Ai  A3(pT'-‘  -  (2.24) 

Indeed,  the  term 


C*=  At^p^-  A3$^$ 


i#-l 


jjrj 


(2.25) 


can  be  thought  of  as  a  conservation  correction  to  the  equation.  This  term  affects  both  the  mass  conservation 
and  time  accuracy  of  the  equation.  Substituting  equation  (2.21)  into  equation  (2.20)  now  yields 


Jji ) 


N 


($^+i  -  $J")  +  ($J"  -  ^ 


iiVx 


.N  ^N-U 


jN 


.#-1-1 


i# 


rN, 


.#-1-1 


i# 


iV+1 


IP  fp 


„#-l  i^N-1  ^N-2 


+  C*(p^^p 


2)  =  0 

(2.26) 


The  density  in  the  space  terms  in  equation  (2.26)  can  be  computed  at  the  N  time  step  with  an  error  of  only 
At.  For  convenience  the  space  terms  are  written  in  “delta”  form.  For  example,  the  streamwise  flux  term  is 
written 

/  „  \  #-|-l  //i\#  /  n  \  ^ 

^afu)  =SJ{-)  Ai5,($^+l-$^)  +  5jfu 


y 


(2.27) 


After  applying  the  delta  form,  dividing  through  by  — (/dAt^),  and  collecting  terms,  equation  (2.26)  becomes 


I  +  At^Ai$f5.+ At^^A3$„5„-AT^AT^  +  ‘  ( 


-  At""  At 

J 


#  Ar#+1  /  ” 


,# 


#Ae#+l  (  ” 


I)  „#; 


0^ 


J^^[j  I 


#A 


#ArAr-hl 


+ 


V  , 

N-1' 


At""  At 


.1  \  .J.V 


At 


# 


—  ($^  _  2$""“"  -|-  $""“"")  -|-  ($""  -  $""“"-) 


,#-l  ,  ^#-2 


#  ^#-A 


+ 


P 


P 


N 


\At^-^ 
At^ 


Ai$f-l5^+ A3$^-l5^ 


($#  _  $#-l)  +  (  AA  )  At"(p"  -  p"-i)  (2.28) 


N(-N  '^N-l\ 


Equation  (2.28)  may  be  written  in  the  following  compact  form: 

L($^+l  -  $^)=  R*(p^,$^)-H($^-  $^“^)  +  C*(p^,p"' 


$  V 


,/?",/?"  “^)  (2.29) 
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whe 


and 


J 


R*  =  (  ^  ) 


N 


(2.30) 


°‘=  (^)  (^)(*'''-2*'''-‘+f''-^)+(|^)A<'''(.''-.'''-‘, 


/'/3"“'\  . 


N 


(2.31) 


The  left-hand  side  term  L  is  simply  the  time  update  term  which  arises  from  the  delta  form  of  the  equation.  The 
term  A$^  actually  arises  from  the  left-hand  side  operator  but  is  placed  on  the  right-hand  side  for  operational 
reasons ,  R*  is  the  spatial  flux  evaluated  at  time  step  N,  and  C*  is  the  conservation  correction  term. 


2.6.  Implementing  the  Algorithm 

2.6.1.  Basic  Difference  Operators 


Equation  (2.28)  must  now  be  implemented.  The  basic  difference  operator  for  time  is 


$iV+l  _ 
At^+1 


= 


(2.32) 


The  velocity  term  which  is  used  in  the  computation  of  L  and  R,  is  formed  with  a  central  difference 


The  computation  of  the  corresponding  term  must  allow  for  jumps  in  $  across  branch  cuts.  For  convenience, 
the  grid  is  assumed  to  be  permeated  with  horizontal  branch  cuts  which  lie  slightly  above  the  grid  points.  Each 
grid  point  has  a  jump  in  potential  E  associated  with  it.  The  value  of  F  is  zero  everywhere  except  along  an 
actual  branch  cut.  A  difference  expression  which  accounts  for  this  field  of  branch  cuts  is 

($,),  =  -  ($,  +  r^.)  +  -  ($^._i  +  r^._i)]  (2.34) 


Figure  5  illustrates  the  velocity  computation  near  a  vortex  branch  cut. 

2.6.2.  Flux  Operator  Differencing 

Each  of  the  spatial  flux  terms  is  made  up  of  the  product  of  three  terms:  (1)  a  velocity  term,  (2)  a  density 
term,  and  (3)  a  metric  term.  This  product,  which  is  the  local  mass  flux,  is  computed  at  the  midpoint  between 
two  grid  nodes.  The  velocity  term  is  computed  by  using  a  one-sided  difference 

($5)i+(l/2)  =  $i+l-$i  (2.35) 


2.6. 2.1.  Flux  density  term  (switching) .  The  density  terms  are  used  to  aid  the  stability  of  the  algorithm. 
Because  the  HSI  problem  is  transonic,  the  type  of  the  equation  will  change  from  elliptic  to  hyperbolic  depending 
on  local  Mach  number.  Stability  is  achieved  by  switching  the  type  of  operator  with  the  equation  type.  To 
illustrate  the  requirement  for  switching,  consider  the  following  two-dimensional  equation: 

^  XX  ^  yy  =  0  (2.36) 
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Figure  5.  Computation  of  velocity  at  grid  point  in  field  permeated  with  branch  cuts. 


We  define  operators  such  that 


-  ^i-1 
Ax 


(2.37a) 


(2.37b) 

(2.38) 


is  stable  for  Moo  <  1,  and  that 

(l-M^)v^V^$+Vj,Aj,$  =  0  (2.39) 

is  stable  for  Moo  >  1.  A  problem  arises  when  we  try  to  use  either  of  these  schemes  alone  for  a  transonic  flow. 
Equation  (2.39)  is  divergent  if  Moo  <  1  and  equation  (2.38)  is  convergent  for  Moo  >  1  only  if 

Aj. 

(l-M^)A, 

which  is  impractical  for  Moo  near  1.  Stability  is  achieved  by  switching  from  equation  (2.38)  to  equation  (2.39) 
depending  on  local  Mach  number.  In  the  present  algorithm,  this  switching  would  lead  to  a  complex  matrix 
form  which  is  costly  to  evaluate.  However,  Holst  and  Ballhaus  (ref.  19)  introduced  another  type  of  switching 
which  is  well  suited  to  the  conservative  form  of  the  equation  and  reduces  the  complexity  of  the  matrix.  In 
their  scheme,  the  density  is  evaluated  centered  in  midcell  for  subsonic  flow.  For  supersonic  flow  the  density  is 
upwind  biased,  that  is,  evaluated  at  the  upstream  cell  center.  (See  sketch  A.)  A  parameter  is  then  employed 
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Spacial  differencing  and  stablility: 


Stable  difference  for  spacial  flux  in  subsonic  flow  is 

(p^P^  =  V(Pi^l/2A(^)  ^ 


<Pi+l/2=7(Pi  +  Pi+l) 


Centered  difference: 

Stability  in  supersonic  regions  requires  upwing  biasing  and 
a  simple  way  to  accomplish  this  is 


(p^^)^  =  V(Pi^l/2A(^)  ^  Cj 


Pi-3/2  Pi-1/2 

This  remains  a  centered  scheme;  p  is  shifted  however 

Sketch  A 


which  switches  the  density  from  centered  to  upstream  based  on  the  local  Mach  number.  The  flux  (without  the 
metric  term  for  simplicity)  is 


(1  -  Vi) - - -  +  Vi - - - 


n  Pi-\  ,  Pi-\^  Pi-‘l 

(1  -  ^ - +  ^i-1  - ^ - 


(^i+1  -  ^i) 

($i-  $i-l) 


(2.40) 


The  switching  parameter  v  is  defined  as 


o 

1 

II 

(1<C  <  10) 

v=  0 

(v*  <  0  (subsonic)) 

V  =  1 

(v*  >  1  (supersonic)) 

(2.41) 


The  parameter  C  is  used  to  provide  additional  numerical  viscosity  in  the  supersonic  region.  The  form  of  C  is 
completely  arbitrary.  In  the  present  method,  C  is  varied  linearly  with  the  local  Mach  number  by  using 
the  following  equation: 

C  =  WMi  -  10.8 

where  this  form  was  determined  by  numerical  experiment.  The  critical  density  p*  is  determined  with  the 
free-stream  Mach  number  as 

2  A  j-1^  Nll/(T-1) 


P*  = 


Lt  + 1 


1  + 


-M. 


2.6. 2.2.  Flux  metric  term  (con, si, stent  differencing).  Computation  of  the  flux  metric  terms  poses  a  special 
problem.  The  flux  metric  term  is 


J 


n  /  i+  1/2 


_  (e+3i 

Jr 


(2.42) 


i+ 1/2 


11 


The  simplest  method  of  evaluating  this  is  to  average  the  values  at  the  nodes  so  that 


i+1/2 


(2.43) 


This  evaluation  can  lead  to  an  error  for  a  grid  with  stretching.  The  error  occurs  because  the  values  and  rfy 
are  computed  by  using  central  differences  (eqs.  (2.13)  and  (2.14)),  and  therefore,  information  from  node  {i  —  1) 
to  (i  +  2)  is  used  to  compute  the  metric  term  at  (f  +  1/2).  This  “extra”  data  acts  to  reduce  the  accuracy  of 
the  metric  computation  especially  for  regions  where  the  metrics  change  rapidly.  The  effect  of  the  diminished 
accuracy  is  to  introduce  an  artihcial  mass  into  the  flux  computation.  The  extra  mass  can  be  removed  by  the 
use  of  a  free-stream  subtraction  matrix.  This  matrix  is  generated  by  specifying  free-stream  conditions  on  the 
mesh  ($  =  M^x)  and  computing  the  value  of  R^o  =  R($  =  Moo  x).  The  matrix  R^^  is  then  subtracted  from 
R  in  subsequent  computations  to  restore  the  mass  balance. 

Flores  et  al.  (ref.  20)  recently  proposed  a  superior  method  of  computing  the  metric  which  eliminates  this 
problem.  Their  method  involves  calculating  the  flux  metric  terms  at  the  same  points  in  space  at  which  the 
flux  is  differenced.  The  method  is  referred  to  as  the  “consistent-metric  method.”  Consider,  for  example,  the 
incompressible  mass  conservation  equation  which  must  hold  for  the  free-stream  subtraction  to  be  zero 


Here,  the  density  is  set  to  1,  and  from  the  chain  rule. 


=  ^($X  -  Vx^ri)  =  -  Vy^ri) 


$  — 
*  71  - 


n  =  — (^x  -  e^^j)  =  -  ^y^i) 

'lx  <^y 


(2.44) 


(2.45) 


With  equations  (2.45)  and  (2.13),  equation  (2.44)  becomes,  for  free-stream  conditions. 


MooiVy^-  )  =  0 


(2.46) 


It  therefore  follows  that  the  metric  difference  operators  must  commute  in  order  to  produce  a  zero  free-stream 
subtraction.  For  this  method,  the  primitive  metrics  (in  the  /  direction)  are  computed  by  using  a  backward 
difference 


the  derivatives  are 


~  yi+1  yi 


y^i+iiij 


1 

2 

1 

2 


/ xj_^l  a;  ■_!  \ 

V  2 

/ t/j+1  -  yj-\  \ 

V  2  A-+1 


■ 


The  determinant  Dm  is 


Dm  —  x^yrj  —  X 


(2.47) 


(2.48) 


(2.49) 
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and 


The  flux  metric  term  becomes 


ix 


f]y  = 


J  m 


D 


-1 

m 


aa  _  il  +  il 

Jn/i+l/2 


(2.50) 


(2.51) 


The  metric  term  (  4^ )  is  computed  in  an  analogous  fashion. 

VTr/j+i/2 

The  consistent  difference  therefore  localizes  the  computation  of  the  metric;  thus,  the  extraneous  data  are 
eliminated  from  the  computation.  Therefore,  three  sets  of  metrics  exist:  one  each  for  the  ^  and  rj  flux  terms 
and  one  at  the  nodes.  Although  the  consistent  metric  method  requires  more  storage  space  and  complicates  the 
coding,  it  provides  a  more  robust  method  of  computing  the  metrics  and  has  been  incorporated  in  the  present 
method.  The  final  form  of  the  flux  term  is 


M  ^  ^2+  1  +  Pi  ,  Pi  +  Pi- 1 
(1  -  ^i) - 5 - +  - 5 - 


Ai  , 

($i+l-$i) 


n/ i  +  1/2 


n  ,,  \Pi+Pi-l  I  ,,  Pi-l  +  Pi-2 
{  ^  -  ^i-l) - - r  - r - 


(2.52) 


n  J  i-1  /2 


Equations  (2.52),  (2.41),  (2.40),  (2.35),  (2.34),  (2.33),  and  (2.32)  are  used  in  the  implementation  of 
equation  (2.28). 


2.6.3.  Approximate  Factorization 

A  noncompact  matrix  inversion  as  follows  is  still  required  to  solve  equation  (2.29): 

LA$^+1  = 


In  order  to  reduce  the  matrix  to  compact  form,  the  operator  L  is  approximately  factored.  The  equation 
becomes 

(2.53) 

The  operators  and  are  chosen  so  that  (1)  their  product  is  approximately  equal  to  L  (to  within  an  error 
which  does  not  exceed  the  discretization  error),  (2)  only  simple  matrix  operations  are  required  to  obtain  the 
solution,  and  (3)  the  overall  scheme  is  stable.  The  present  method  uses  factors  which  are  associated  with  the 
two  space  derivatives,  and  this  leads  to  an  ADI  type  scheme.  The  factor  is 


1+  Af^A3$^(5^  -At^At^^^ 


J  n 

w 


As 


N  n 

P  dn 


(2.54) 


13 


The  form  of  is  similar.  The  final  form  of  equation  (2.29)  is 


I  +  -  A(«A(«+1  (^)  0„  (^)  p-^d.. 


j3N 


Af^Af^+1 

+  ($^-  + 

f  pN-V 


b^  (^)  +  in  (^)  P^in^^ 


j3N 


+ 


/3^ 


At 


# 


Ai$f-l6j+ A3$f-% 


($iV  _  ^#-1^  ^  _  ^#-1^  ^2.55) 


Both  operators  L^  and  Ljj  yield  tridiagonal  matrices.  For  example, 


=  A$*_i  + 


(2.56) 


whe 


A  = 


-At^ -  At^At^+^  ^  — 


f]N 


A3\  Pj  +  Pj-\ 

Jji  J  j-\l‘2  2 


1  +  At^At^+^  -7^  1  (  ^ 


N 


+At^At^+^  ( ^ 

\P^ 


A3 

J 


a.7'+l  +aj 

2 


Jji  /  j  +  1/2  2 


C  = 


Ji/ j-1/2 

At^A3$^  _  ^ 

2 


A3 

J  n 


Pj+\  + 


i+1/2 


(2.57) 


2.6.4-  Implementing  Boundary  Condition, s 


The  boundary  conditions  are  implemented  implicitly  in  the  algorithm.  This  implemention  will  require 
modifications  to  both  the  right-  and  left-hand  sides  of  the  rj  sweep  of  the  equation.  For  the  upper  surface  of 
the  airfoil,  the  operator  becomes 


At^+lA3<F^A5„  At^At^+l 


P^  [\JnJj+il/2)  2 


n  /  S 


(2.58) 


Here  the  term  A5'u  refers  to  the  change  in  velocity  on  the  upper  surface  from  time  A  to  A  +  1.  Because  this 
is  a  known  quantity,  it  can  be  brought  to  the  right-hand  side.  The  equation  then  becomes 


‘^j+l 


iV-1  AliV  .  ^^#-1 


A$J_i  +B$*+  =  R+  A5^"“‘AA"A3$ 


n 


-  At^As^.AS:  -  At-  At-+‘  (  ^  M  )  PS  AS: 


.N  A  ,#  Aitv-hl  / 


A3 


•n  /  S 


(2.59) 
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where 


A  =  0 


B  = 


At^+^At^  fA^\  Pj+^  +  pj 

VjJi+i/2  2 


Q  -  _  ^j+1  +  Pi 

Wn/j+l/2  2 

R  represents  that  portion  of  the  right-hand  side  which  is  unchanged  at  the  boundary,  and  the  subscript  u  refers 
to  the  upper  surface  conditions.  The  modification  on  the  lower  surface  is  similar  to  that  on  the  upper  surface 
and  can  be  determined  by  symmetry. 


2.6.5.  Solution  Steps 

The  solution  is  obtained  in  three  steps  by  the  following  equations: 


=  R 

(2.60) 

(2.61) 

$  tV-l- 1  =  $  ^  q.  1 

(2.62) 

The  value  is  the  solution  of  the  full- potential  equation  for  the  flow  about  the  airfoil  at  the  next 

time  step.  Given  this  solution,  the  velocity  and  pressure  on  the  airfoil  surface  may  be  calculated.  The  vortex 
model  is  introduced  as  a  modification  to  either  the  boundary  conditions  or  the  basic  algorithm. 


2.7.  Steady-State  Algorithm 


A  special  form  of  the  algorithm  is  employed  in  the  solution  of  a  steady  problem.  For  steady  problems, 
equation  (2.55)  is  modified  to  remove  all  temporal  terms  and  the  resulting  equation  is  solved  with  pseudo  time 
terms  which  act  to  update  the  solution.  Equation  (2.55)  becomes 


1-At^ 


drj 


N 


(2.63) 


The  value  of  At^  is  then  oscillated  for  a  number  of  time  steps.  Each  successive  value  of  At  acts  to  reduce  the 
magnitude  of  the  error  in  a  limited  frequency  range.  By  oscillating  At,  the  error  for  a  wide  range  of  frequencies 
is  reduced  to  make  the  most  efficient  use  of  each  computational  sweep.  With  this  method,  approximately  400 
“time  steps”  are  required  to  drive  the  residual  to  an  acceptable  value,  whereas  2000  steps  are  required  with 
the  full  algorithm  (with  all  the  time  terms  included). 

The  residual  is  determined  by  the  following  steps: 

1.  Compute  the  value  of  R  at  the  first  step 

2.  Survey  R  to  obtain  its  maximum  local  value  R| 

3.  At  each  subsequent  computational  step,  obtain  R„  (the  maximum  local  value  of  R) 

4.  Determine  the  value  of  Rn/R^  which  is  the  normalized  residual 
When  this  value  reaches  10”"^  the  computation  has  converged. 
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3.  Vortex  Models 

In  section  2,  an  algorithm  for  solving  for  the 
potential  flow  field  around  an  airfoil  in  transonic 
flow  was  presented.  This  algorithm  is  now  modified 
to  include  the  effect  of  a  two-dimensional  vortex 
passing  near  the  airfoil.  The  four  models  for  the 
vortex  discussed  in  section  1.3  are  used:  the  angle- 
of-attack  method,  the  lifting-surface  method,  the 
branch-cut  method,  and  the  split-potential  method. 
These  methods  may  be  grouped  into  two  categories. 
(See  fig.  6.) 

The  first  category  called  the  surface-specification 
methods  model  the  effect  of  the  vortex  as  an  imposed 
normal  velocity  distribution  on  the  airfoil  surface. 
Both  the  angle-of-attack  and  lifting-surface  methods 
fall  into  this  category.  These  methods  originated 
within  linear-integral  rotor  theories.  The  effect  of 
the  vortex  on  the  general  flow  field  is  usually  not 
considered  in  these  theories.  These  models  are  valid 
for  linear  flow  fields. 

For  problems  characterized  by  the  transonic  non¬ 
linearity  (that  is,  with  a  speed  of  sound  that  varies 


throughout  the  flow  field),  a  surface  effect  cannot 
completely  model  the  effect  of  the  vortex;  therefore, 
it  is  necessary  to  insert  the  vortex  explicitly  into  the 
grid.  This  category  is  called  explicit  models;  both 
the  branch- cut  and  split- potential  models  fall  into 
this  category. 

The  vortex  modeling  begins  with  the  ideal  two- 
dimensional  vortex  potential: 

G=L0  (3.1) 

Ztt 

where  9  is  the  angle  subtended  by  the  vortex  and  the 
field  point.  The  tangential  velocity  at  the  field  point 
is 


The  singularity  at  r  =  0  is  the  source  of  numerical 
instabilities  and  requires  the  use  of  an  artificial  core. 
In  the  following  tasks,  the  model  developed  by  Scully 
(ref.  21)  given  in  the  following  equation  is  used: 


Vortex  potential 
or  velocity  field 
specification 


Explicit  models 

Figure  6.  Principal  vortex  models. 


Blade  inflow 
specification 


Surface-specification  model 
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where  a  is  the  vortex  “core”  radius.  The  word 
core  is  used  in  the  conventional  sense  that  is  used  in 
the  helicopter  industry;  that  is,  it  refers  to  a  region 
of  rotational  flow  near  the  vortex  center  of  rotation. 
This  rotational  flow  may  be  restricted  to  a  discrete 
region  in  some  models  or  it  may  be  modeled  by  a 
decay  function  as  in  equation  (3.3).  The  radius  a 
defines  the  region  in  which  the  flow  is  rotational; 
within  this  region,  the  potential  equation  model  of 
the  vortex  is  invalid. 


the  same  time;  such  a  solution  would  be  a  close 
approximation  to  the  exact  solution  of  the  flow  field. 
However,  the  effect  of  a  finite  signal  propagation 
speed  is  still  violated  by  these  methods  and  this 
becomes  increasingly  important  as  the  vortex  nears 
the  airfoil  particularly  for  transonic  flow  conditions. 
The  surface  specification  models  are  not  capable  of 
modeling  the  “time  lag”  between  a  signal  arriving  at 
a  point  on  the  airfoil  surface  near  the  airfoil  and  the 
signal  arriving  on  the  opposite  side  of  the  airfoil. 


The  vortex  is  moved  through  the  computational 
grid  by  integrating  the  flow  velocity  at  the  vortex 
over  the  current  time  step 


Xf +1  =  Xf  +  17,  1 
y/+i  =  y/  +  ] 


(3.4a) 


The  vortex  convection  velocities  (7,  and  17,  can  be 
determined  by  three  different  methods:  (1)  a  priori 
specification,  (2)  interpolation,  and  (3)  the  velocity 
field  of  a  point  vortex. 

By  far,  the  easiest  method  is  to  specify  an  initial 
position  and  then  allow  the  vortex  to  convect  at  the 
free-stream  speed.  The  equations  most  often  used 
are 

y,  =  Moo  1 

17,  =  0  J 

Equations  (3.4b)  produce  a  “fixed-path”  interaction. 
Although  specifying  the  vortex  is  a  trivial  matter 
in  a  2-D  flow  problem,  it  is  the  usual  procedure 
in  3-D  integral  computations  of  advancing  rotors 
because  of  the  cost  and  complexity  of  finding  the 
wake  deformation.  This  method  is  also  very  useful 
for  various  comparison  studies. 


(3.4b) 


3.1.  The  Surface-Specification  Methods 


3.1.1.  Angle -of- Attack  Method 


The  angle-of-attack  method  is  the  simplest  pos¬ 
sible  model  of  the  effect  of  a  vortex  on  an  airfoil. 
Equation  (3.3)  is  used  to  compute  the  velocity  at  the 
airfoil  quarter-chord.  With  this  velocity,  a  vortex- 
induced  angle  of  attack  is  computed.  (See  fig.  7.) 
The  velocity  perpendicular  to  the  chord  line  is  (if 
the  leading  edge  is  at  *  =  0) 


r  70.25 -A, 

V±=  Vg  cos9  =  —  i  — 2^2 
Ztt  \ 

The  vortex-induced  angle  of  attack  is 

.-1 


=  tan 


Co 


(3.5) 


(3.6) 


This  angle  is  added  to  the  airfoil  angle  of  attack.  The 
potential  field  is  then  computed  as  before.  The  angle 
of  attack  is  updated  at  each  time  step  as  the  vortex 
moves  by  the  airfoil. 

The  airfoil  is  therefore  assumed  to  be  a  point  in 
space.  In  order  for  this  solution  to  be  valid,  the  vor¬ 
tex  must  be  far  enough  away  for  this  approximation 
to  be  accurate  (e.g.,  the  signal  must  arrive  at  every 
point  on  the  airfoil  simultaneously). 


The  surface-specification  models  are  produced  by 
modifying  the  airfoil  surface  boundary  conditions 
based  on  an  assumed  vortex  velocity  field.  This 
approach  is  basically  the  same  as  assuming  that  the 
vortex  velocity  may  be  superimposed  on  the  general 
flow  problem  in  the  same  way  as  a  gust  velocity 
would  be  modeled.  The  accuracy  of  this  assumption 
depends  on  the  location  of  the  vortex  with  respect 
to  the  airfoil.  If  the  vortex  is  far  enough  away,  the 
field  which  it  produces  does  resemble  a  gust  field.  A 
constant  velocity  field  produces  the  angle-of-attack 
method,  and  a  variable  velocity  field  produces  the 
lifting-surface  method.  Furthermore,  if  the  vortex  is 
far  enough  away  from  the  airfoil,  the  signal  arrives  at 
the  various  points  along  the  airfoil  at  approximately 
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3.1.2.  Lifting-Surface  Method 


The  lifting-surface  method  is  the  most  general 
form  of  the  surface-specification  models.  As  with 
the  angle-of-attack  method,  equation  (3.3)  is  used  to 
compute  the  induced  velocity  at  the  airfoil  surface. 
(See  fig.  8.)  However,  unlike  the  angle-of-attack 
method,  the  velocity  is  allowed  to  vary  over  the 
surface.  This  provides  a  more  physically  accurate 
model  of  the  effect.  For  the  lifting- surface  method, 
equation  (3.5)  is  modified  and  is 


Figure  9.  Branch-cut  vortex  model. 


Vg  COS  9  = 


{x/c)  -  Xy 

J.2  _|.  ^2 


(3.7) 


With  this  change,  the  computation  proceeds  as  in  the 
angle-of-attack  method. 

The  airfoil  is  therefore  assumed  to  be  a  lifting 
surface.  In  order  for  this  solution  to  be  accurate, 
the  signal  must  arrive  at  the  upper  and  lower  surface 
simultaneously. 


3.2.  Explicit  Models 


3.2.1.  Branch-Cut  Method 


Caradonna  (ref.  22)  was  the  first  to  use  an  explicit 
method  vortex  model  with  a  finite- difference  rotor 
computation.  The  branch-cut  method,  which  he  used 
for  steady  3-D  flows,  is  based  upon  the  known  poten¬ 
tial  solution  for  a  two-dimensional  vortex  (eq.  (3.1). 
This  potential  is  implemented  by  means  of  a  branch 
cut  which  extends  from  the  center  of  the  vortex  to 
the  aft  face  of  the  computational  grid.  (See  fig.  9.) 
A  jump  in  potential  equal  to  F  is  imposed  across  the 
cut.  Because  equation  (2.34)  has  already  been  im¬ 
plemented  to  account  for  the  wake  cut,  no  changes 
in  the  algorithm  are  required. 


At  first  the  branch-cut  method  seems  to  be  well 
suited  to  a  potential  finite-difference  algorithm.  Dif¬ 
ficulties  arise,  however,  in  unsteady  problems  when¬ 
ever  the  vortex  is  moved.  As  the  edge  of  the  cut 
moves  past  a  node,  an  abrupt  change  in  the  local 
potential  occurs.  This  sharp  change  causes  spurious 
waves  which  affect  the  entire  flow  field.  The  problem 
can  be  solved  by  spreading  the  edge  of  the  branch 
cut,  that  is,  by  distributing  vorticity  on  the  vari¬ 
ous  nodes  which  surround  the  vortex  center.  The 
simplest  distribution  involves  the  use  of  the  near¬ 
est  four  grid  points.  The  distribution  is  weighted 
so  that  the  “center  of  gravity”  of  the  vorticity  repre¬ 
sents  the  center  of  the  vortex.  With  four  grid  points, 
this  will  uniquely  determine  the  vorticity  distribu¬ 
tion.  Increasing  the  number  of  points  would  require 
an  arbitrary  distribution  to  be  imposed  upon  the  vor¬ 
ticity.  With  this  modification,  the  vortex  may  be 
moved  from  cell  to  cell  smoothly  and  the  spurious 
waves  are  reduced  (not  eliminated).  The  method  of 
distributing  the  vorticity  is  illustrated  in  figure  10. 
The  vertical  distribution  of  vorticity  is 

Fi  =  F,  (^1  -  (3.8) 


Figure  8.  Lifting-surface  vortex  model. 


and 

F2  =  F,-Fi  (3.9) 

The  horizontal  distribution  of  vorticity  is 

X  1  X 

F4  =  Fi^^  3.10 

X* 

and 


The  main  vortex  at  (xy,  t/^)  is  then  modeled  with  two 
branch  cuts  of  varying  strength;  this  can  be  called  a 
two-cut  model.  The  vortex  may,  in  fact,  be  modeled 
by  any  arbitrary  distribution.  Stremel  (ref.  23)  uses 
a  method  in  which  the  vortex  is  modeled  with  an 
area-weighted  distribution  of  vorticity.  A  parabolic 
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Figure  10.  Distribution  of  vorticity  to  four  points  nearest 

vortex. 

distribution  of  vorticity  in  the  horizontal  direction 
coupled  with  a  linear  vertical  variation  was  used  in 
reference  24.  The  effect  of  distributing  the  branch 
cuts  is  to  create  an  artificial  core  for  the  vortex.  The 
efficacy  of  the  core  is  dependent  upon  the  distribution 
of  the  nodes  which  are  in  the  vicinity  of  the  vortex. 
Because  each  of  the  separate  branch  cuts  represents 
a  separate  subvortex  and  each  subvortex  has  its  own 
singular  point,  the  core  is  very  sensitive  to  the  grid 
geometry. 

Another  problem  associated  with  moving  the  vor¬ 
tex  is  concerned  with  computing  the  vortex  convec¬ 
tion  velocity.  Interpolation  of  the  local  velocities  near 
the  vortex  is  the  only  available  means  of  comput¬ 
ing  the  vortex  velocity  directly.  The  interpolation  is 
complicated  by  the  fact  that  the  vortex  creates  such 
a  large  local  disturbance;  separating  the  effect  of  the 
flow  field  on  the  vortex  from  the  effect  of  the  vortex 
on  the  flow  field  is  difficult.  One  major  shortcoming 
of  the  branch-cut  method  is  no  good  way  exists  to 
separate  the  effects  because  the  branch  cut  contains 
the  combined  potential  of  the  vortex,  the  free  stream, 
and  the  airfoil.  The  vortex -induced  velocities  domi¬ 
nate  the  flow  near  the  vortex  and  make  an  accurate 
interpolation  very  difficult.  The  velocity  at  the  vor¬ 
tex  may  be  approximated  by  using  the  velocity  field 
of  a  point  vortex  in  conjunction  with  the  lift  on  the 
airfoil.  In  this  approach,  equation  (3.3)  is  essentially 
used  in  “reverse”  with  F  being  the  jump  in  potential 


at  the  airfoil  trailing  edge  (and  hence  a  measure  of 
the  airfoil  lift). 

Another  feature  of  the  branch-cut  method  is  the 
fact  that  it  requires  a  difference  equation  (eq.  (2.34)) 
to  implement  the  effect  of  the  vortex.  That  is,  the 
vortex  effect  is  specified  entirely  by  the  potential 
jump  which  is  represented  by  the  differencing  across 
the  branch  cut.  The  accuracy  of  this  difference  is  also 
dependent  upon  the  local  grid  geometry.  Therefore 
the  accuracy  of  the  vortex  model  changes  as  the 
vortex  moves  through  the  mesh.  The  distribution  of 
the  vorticity  on  the  mesh  further  distorts  the  model 
by  increasing  the  mesh  dependence.  A  successful 
branch-cut  model  is  therefore  a  compromise  between 
an  effective  core  model  and  an  accurate  vortex  model. 

3.2.2.  Split- Potential  Method 

An  alternative  to  the  branch-cut  method  is  the 
split-potential  method.  (See  fig.  11.)  In  this  ap¬ 
proach,  the  velocity  is  assumed  to  be  a  combination 
of  a  known  velocity  and  a  perturbation  velocity  as 
follows: 

q  =  V</)+VG  (3.12) 

where  Vq  is  the  known  velocity  field  solution  and  V (j) 
is  a  perturbation  velocity,  which  need  not  be  small. 
The  total  potential  for  the  blade- vortex  problem  can 
be  split  between  the  perturbation  potential  (associ¬ 
ated  with  the  airfoil)  and  the  potential  G,  which  de¬ 
scribes  the  vortex  velocity  field,  as  follows: 

^=(f)+G  (3.13) 

Any  potential  algorithm  may  be  modified  in  this  way 
to  include  the  effects  of  a  known  velocity  component 
and  a  perturbation  velocity.  Furthermore,  the  po¬ 
tential  G  need  not  represent  a  vortex  but  can  in  fact 
represent  any  flow  field  which  independently  satis¬ 
fies  the  potential  equation.  When  equation  (3.13) 


Figure  11.  Split- potential  vortex  model. 
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is  applied  to  equation  (2.28),  the  following  equation  Ames  Research  Center.  The  first  step  is  to  recast 

results  (in  the  compact  form  similar  to  eq.  (2.29)):  equation  (3.14)  as 


L(/,  =  R*(/', 

+  c  *(p^  4>^-^ , 

+  C*(G^ ,G^-\  ^ ^ pN- 1) 

+  A<?i"  +  AG"-L(p"  AG^+1)  (3.14) 


L(p",  A$"+‘)  =R.*(p- 


N+l  \  _  ■ 


„N  iN 


+  G") 


,  r^*/  N  N-l  ,N  ,N-1  iN-i 

+  C  (fi  ,p  ,<t>  ,<t>  ,<t>  ,13  ,13  ) 

+  C*(G",  G^-\  G^-‘^,/3^,  /3"-i) 

+  A<?i"  +  AG"  (3.18) 


The  Bernoulli  equation  undergoes  a  similar 
modification 

+  G^)  (3.15) 

The  left-hand  side  of  equation  (3.14)  is  identical 
to  the  original  algorithm  (eq.  (2.28)).  The  right- 
hand  side  contains  additional  spatial  and  temporal 
gradient  terms  in  G,  including  the  update  operator 

L(G). 

Implementation  of  equation  (3.14)  has  proven 
to  be  a  challenge  to  several  researchers  who  have 
sought  to  simplify  the  equation.  (See  refs.  7,  13, 
and  17.)  The  principal  focus  of  these  studies  has 
been  in  eliminating  the  temporal  gradient  terms  in 
G.  These  terms  pose  a  particular  problem  because 
they  involve  the  potential  G  explicitly.  Computing 
these  terms  requires  the  tracking  of  a  branch  cut 
through  the  flow,  in  effect,  the  split-potential  model 
is  reduced  to  a  branch-cut  model.  This  will  be 
particularly  difficult  for  the  complex  geometry  of  the 
full  3-D  problem.  Furthermore,  it  is  advantageous  to 
minimize  the  computational  requirements  as  much 
as  possible.  McCroskey  and  Goorjian  (ref.  13)  have 
shown  (for  a  small  disturbance  formulation)  that 
L(G)  and  AG^  can  be  eliminated  since  the  vortex 
potential  (eq.  (3.1))  is  a  solution  to 

LAG  =  0  (3.16) 

Therefore,  a  “small  disturbance  version”  of  equa¬ 
tion  (3.14)  would  be 

L{p^ ,  =  R*(/9^,  +  G^)  +  A(p^  (3.17) 

because  C*  terms  are  not  present  in  a  small  distur¬ 
bance  form.  Sankar  and  Malone  (ref.  17)  restricted 
the  solution  to  his  algorithm  to  a  so-called  “weak 
split-potential”  approach  in  which  the  temporal  gra¬ 
dient  and  most  of  the  spatial  terms  in  G  are  simply 
dropped. 

The  present  method  is  neither  a  small  disturbance 
or  weak  split-potential  form.  In  spite  of  this,  the 
algorithm  can  be  simplified  by  using  a  method  pro¬ 
posed  by  Roger  Strawn,  Aeroflightdynamics  Direc¬ 
torate,  U.S.  Army  Aviation  and  Missile  Command, 


The  update  operators  L  have  been  recombined  to 
include  the  total  potential.  The  temporal  conser¬ 
vation  correction  term  C*  in  G  has  been  retained, 
since  there  is  no  way  to  effectively  separate  the  G 
and  (j)  parts  of  the  density.  The  effect  of  these 
terms  on  the  solution  is  discussed  in  section  4.  So¬ 
lution  of  equation  (3.18)  is  achieved  as  before  with 
equations  (2.60),  (2.61),  and  (2.62),  which  produces 
^A-l-C  A  final  step  is  added  to  obtain  (j) 

^  _  ag^+^  (3.19) 

In  implementing  these  equations,  the  velocity  due 
to  the  vortex  is  computed  by  equation  (3.3).  The 
velocity  components  are 


1 

II 

0 

f  y-  yv\ 

\r2  +  / 

(3.20) 

G 

^  X  —  Xy\ 

(3.21) 

27r\ 

+  a?) 

These  equations  are  transformed  into  G^  and  G^  by 
G^  =  GxX^  +  GyUf:  (3.22) 

and 

Gj]  =  GxXrj+  GyUfj  (3.23) 

The  time  derivative  terms  are  obtained  with  the 
chain  rule  as 

Gx  =  G^^x  +  Gfjrjx  +  Gttx  (3.24) 

This  value  may  be  determined  most  easily  in  an  axis 
system  fixed  to  the  vortex.  In  this  system,  the  vortex 
is  fixed  and  the  airfoil,  to  which  the  grid  is  attached, 
moves  past  it.  In  this  system 

G;  =  0  ) 

t]t=  -Vv  \  (3.25) 

=  —Uv  J 

therefore, 

AG  =  At^Gx  (3.26) 
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Use  of  equation  (3.26)  eliminates  the  need  to  com¬ 
pute  G  explicitly  and  hence  the  need  to  track  a 
branch  cut.  Equations  (3.20)  through  (3.26)  are 
used  to  implement  the  split-potential  method  in  the 
present  algorithm. 

4.  Computational  Results 

The  presentation  of  computed  results  is  divided 
into  four  sections.  Section  4.2  reports  results  ob¬ 
tained  for  simple  flows  (with  no  vortex  interaction) 
in  order  to  establish  the  validity  of  the  basic  algo¬ 
rithm.  Section  4. 3  is  an  introduction  to  the  results 
obtained  in  a  HSI  computation.  The  major  features 
of  the  interaction  are  presented  for  a  typical  case. 
The  graphical  presentation  of  data  which  describes 
the  interaction  is  presented  and  explained.  The  sen¬ 
sitivity  of  the  algorithm  to  both  time  step  and  ran¬ 
dom  disturbances  is  discussed.  Section  4.4  presents  a 
comparison  of  the  results  obtained  with  the  four  vor¬ 
tex  models.  Comparisons  with  experimental  data  are 
also  made.  Section  4.5  presents  a  parametric  study 
of  the  effects  of  some  key  parameters  such  as  Mach 
number,  vortex  strength,  miss  distance,  vortex  core 
size,  and  angle  of  attack. 

4.1.  Executing  the  Algorithm 

4.1.1.  Executing  Baste  Airfoil  Solution 


Figure  12.  Maximum  residual  convergence  history  for  NACA 
0012  airfoil  at  several  Mach  numbers  and  a  =  0°. 


A  basic  airfoil  solution  consists  of  the  steady  flow 
around  an  airfoil  at  a  fixed  angle  of  attack.  The 
steady-state  version  of  the  algorithm  (eq.  (2.59))  is 
used  to  generate  this  solution.  Figure  12  depicts 
the  maximum  residual  convergence  history  for  an 
NACA  0012  airfoil  at  a  =  0°  and  several  Mach 
numbers.  The  residual  is  reduced  by  oscillating  the 
pseudo  time  At  between  the  values  0.001  and  5  for 
the  first  20  iteration  steps.  This  oscillation  drops  the 
value  of  Rmax  approximately  1  order  of  magnitude. 
The  pseudo  time  is  then  fixed  at  At  =  0.001  to 
minimize  shock  oscillation  problems  for  the  higher 
Mach  numbers.  Fixing  At  results  in  the  unusual 
flattening  of  the  curve  at  about  100  time  steps.  After 
this  point.  At  =  0.001  is  not  the  optimum  choice 
to  reduce  R.  However,  the  current  scheme  provides 
a  suitable  method  for  reducing  R  to  an  acceptable 
value  over  the  widest  range  of  flow  conditions. 

4.1.2.  Executing  HSI  Problem 

The  solution  of  a  blade-vortex  interaction  prob¬ 
lem  proceeds  in  two  steps.  The  first  step  is  the  com¬ 
putation  of  a  steady-state  solution.  The  steady-state 


solution  may  be  produced  with  or  without  the  pres¬ 
ence  of  the  vortex;  however,  for  convenience  the  vor¬ 
tex  is  included.  The  vortex  is  initially  fixed  at  some 
distance  upstream  (usually  about  10  chords)  and  the 
steady  flow  about  the  airfoil  is  computed.  After  the 
steady-state  solution  is  obtained,  the  vortex  is  al¬ 
lowed  to  move  along  a  path  determined  by  equa¬ 
tions  (3.4a).  A  specified  path  is  produced  by  using 
equations  (3.4b)  to  compute  vortex  velocities.  As  the 
vortex  moves  through  the  computational  grid,  the 
time  step  is  varied  by  using  the  following  relationship 
which  was  established  by  numerical  experimentation: 


At  = 


0.5  cos[7r(*^  +  10.405)] 
10.405 


+  0.522 


Equation  (4.1)  is  designed  to  provide  a  sufficient 
number  of  steps  within  each  grid  cell  to  assure  proper 
resolution  of  the  vortex.  Thus,  when  the  vortex  is  far 
from  the  airfoil.  At  is  large  to  minimize  the  number  of 
steps  required  for  the  solution;  and  when  it  is  near  the 
airfoil.  At  is  small  to  enhance  accuracy.  Application 
of  a  varying  time  step  reduces  the  time  requirements 
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for  computation  by  as  much  as  a  factor  of  6  over  a 
constant  time  step. 


Jf.1.3.  Data  Recovery 

The  velocities  in  the  flow  field  are  computed  by 
the  central  difference  equations  (2.33)  and  (2.34).  In 
the  physical  frame  they  become 


«  =  +  ^yVx 

V  =  ^  r]f]y 


(4.2) 


Pressure  is  calculated  by  using  the  isentropic  flow 
relations  which  lead  to 


Distance  traveled,  chords 


Poo 


(l/2)pU 


(pT  _  1) 


(4.3) 


Lift  on  the  airfoil  is  calculated  by  using  the  trape¬ 
zoid  rule  to  integrate  the  pressures  to  get 


Cl 


2  C  Acp  dC 

Jo 


(4.4) 


where  (  =  ^ x /c.  The  integration  with  respect  to 
the  square  root  of  the  surface  coordinate  increases 
the  accuracy  in  the  leading-edge  region.  The  Mach 
number  is 


a 


2 

oo 


(4.5) 


The  density  is  computed  directly  from  the  Bernoulli 
equation  (eq.  (2.2)) 


P  —  Poo 


^2 


1 


4.2.  Results  From  Basic  Airfoil  Solution 

Before  presenting  the  results  for  the  HSI  problem, 
it  is  useful  to  demonstrate  the  validity  of  the  basic 
algorithm  by  comparisons  with  the  results  obtained 
by  other  researchers  and  by  experiment  for  steady 
and  unsteady  flow  problems. 

4.2.1.  Comparison  With  Original  Algorithm 

The  basic  scheme  used  in  the  present  method 
is  an  extension  of  a  method  developed  by  Steger 
and  Caradonna  (ref.  2).  The  first  step  in  verifying 
the  present  method  is  to  demonstrate  the  ability 
to  reproduce  their  results.  The  predicted  midchord 
values  of  Cp  of  a  parabolic  arc  airfoil  as  it  thickens 
then  thins  for  a  free-stream  Mach  number  of  0.85 


Figure  13.  Pressure  coefEcieuts  at  midchord  for  parabolic 
arc  airfoil  as  it  thickens  then  thins  for  =  0.85  and 
T  =  0.104.  Small  disturbance  grid  used. 


and  maximum  thickness  r  of  0.104  are  presented  in 
figure  13.  The  results  for  the  present  method  were 
predicted  by  using  a  “small- disturbance”  type  mesh 
(see  step  1  in  section  2.2);  reference  2  also  uses  such  a 
mesh.  Figure  13  demonstrates  very  close  agreement 
between  the  present  method  and  reference  2. 

4.2.2.  Subcritical  Flow 

Although  HSI  is,  in  general,  a  transonic  flow  prob¬ 
lem,  the  demonstrated  ability  to  compute  subcritical 
flows  is  necessary  to  establish  confidence  in  the  al¬ 
gorithm.  Figure  14  presents  computed  results  for  an 
NACA  0012  airfoil  at  =  0.63  and  a  =  2°  along 
with  the  results  of  Holst  (ref.  25). 

Figure  15  presents  calculated  results  (performed 
by  Michel  Costes  of  the  Office  National  d’Etudes  et 
de  Recherches  Aerospatiale  (ONERA),  who  had  been 
provided  with  the  present  code)  for  an  Aerospatiale 
RA16SC1  airfoil  at  =  0.30  and  a  =  0°.  The  re¬ 
sults  are  compared  with  experimental  data  generated 
by  ONERA  and  provided  for  this  comparison.  Both 
comparisons  demonstrate  the  ability  of  the  present 
algorithm  to  predict  accurately  the  subcritical  flow 
about  an  airfoil.  The  correlation  with  the  Aerospa¬ 
tiale  RA16SC1  airfoil  is  particularly  good  considering 
the  large  difference  between  the  surface  grid  line  (ob¬ 
tained  from  the  Joukowski  airfoil)  used  to  represent 
the  airfoil  surface  and  the  actual  airfoil  surface. 

4.2.3.  Supercritical  Case 

A  demonstrated  ability  to  compute  supercritical 
flow  is  necessary  in  order  to  proceed  with  the  HSI 
problem.  Figure  16  presents  computed  results  for 
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Figure  14.  Pressure  coefficients  for  NACA  0012  airfoil  at  Moo  =  0.63  and  a  =  2°. 


Figure  15.  Pressure  coefficients  for  Aerospatiale  RA16SC1 
airfoil  at  Afoo  =  0.30  and  a  =  0°.  Data  from  Michel 
Costes,  ONERA. 

an  NACA  0012  airfoil  at  Moo  =  0.8  and  a  =  0°. 
The  surface  pressures  are  compared  with  results 
generated  by  Bridgeman,  Steger,  and  Caradonna 
(ref.  10),  which  were  provided  by  him  for  this  com¬ 
parison.  Close  examination  of  the  figure  reveals  a 


slight  difference  in  shock  location  between  the  two 
models.  This  difference  is  probably  caused  by  the  dif¬ 
ference  in  grid  metric  computation.  The  Bridgeman 
method  uses  the  free-stream  subtraction  method  to 
remove  the  metric  truncation  error  and  the  present 
method  uses  consistent  metric  differencing. 

Figure  17  presents  calculated  results  for  an  Aero¬ 
spatiale  RA16SC1  airfoil  at  =  0.76  and  a  =  0° 
(provided  by  Costes).  The  results  are  compared 
with  experimental  data  generated  by  ONERA.  The 
correlation  is  good,  except  for  a  slight  difference 
in  shock  location  which  is  probably  due  to  viscous 
effects  on  the  real  airfoil. 

4.2.4-  Oscillating  Airfoil 

The  ability  to  predict  unsteady  flows  is  necessary 
in  order  to  proceed  with  the  HSI  problem.  Figure  18 
presents  calculated  results  for  an  NACA  0012  air¬ 
foil  at  Mqo  =  0.755  and  a  =  0°  and  experimental 
results  from  Goorjian  and  Guruswamy  (ref.  26). 
Comparison  of  the  two  results  shows  excellent  agree¬ 
ment.  Figure  19  presents  the  results  for  an  oscil¬ 
lating  NACA  0012  airfoil.  Figure  18  is  the  steady 
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Figure  16.  Pressure  coefficients  for  NACA  0012  airfoil  at 
Afoo  =  0.80  and  a  =  0°. 


Figure  18.  Pressure  coefficients  for  NACA  0012  airfoil  at 
Afoo  =  0.755  and  a  =  0° . 


Exp.  Calc.  Sur. 
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Figure  17.  Pressure  coefficients  for  Aerospatiale  RA16SC1 
airfoil  at  Afoo  =  0.76  and  a  =  0°. 

starting  point.  The  airfoil  is  oscillating  at  a  reduced 
frequency  of  fc  =  ujc/Uqq  =  0.1628  about  a  =  0° 
with  ttmax  =  ±2.51.  The  predicted  pressures  are 
slightly  different  from  the  measured  values  because 
of  the  difficulty  in  exactly  matching  the  test  points. 
These  results  indicate  that  the  algorithm  can  predict 
unsteady  transonic  flow. 

Figure  20  presents  calculated  results  for  an 
Aerospatiale  RA16SC1  airfoil  with  an  oscillating 
25-percent  flap  at  =  0.30  and  a  =  0°.  The  re¬ 
sults  are  compared  with  experimental  data  for  three 


periods  of  motion.  The  pressures  are  presented  in  a 
format  based  on  the  following  equation: 

(4.6) 

k 

where  (c|,)o  is  the  steady  pressure  and  Lp  is  the  phase 
angle  of  the  oscillation.  The  correlation  is  good 
except  for  the  flap  pressures  themselves.  However, 
the  grid  is  sparse  in  this  region  and  the  difference 
could  easily  be  due  to  poor  resolution. 

4.3.  Introduction  to  HSI 

Jf.3.1.  Generic  Vortex  Interaction 

Before  presenting  the  vortex  model  comparison, 
a  discussion  of  the  basic  features  of  the  blade-vortex 
interaction  will  be  useful.  Figure  21  is  a  representa¬ 
tive  plot  of  lift  coefficient  versus  vortex  location  for 
an  NACA  0012  airfoil  at  a  =  0°  and  =  0.60,  and 
a  vortex  with  an  equivalent  lift  of  ^  =  0.400,  and  a 
vertical  miss  distance  of  0.251  chord  moving  along  a 
fixed  path  at  constant  speed.  The  computation  uses 
the  split-potential  method  to  model  the  vortex.  The 
several  points  of  interest  which  occur  during  the  in¬ 
teraction  are  labeled  in  the  figure.  Point  1  is  the  ini¬ 
tial  condition  of  the  airfoil  in  relatively  undistorted 
flow  (the  vortex  is  far  upstream).  The  sign  of  the 
vortex  is  such  that  it  generates  a  downwash  on  the 
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Figure  19.  Pressure  coefficients  for  oscillating  NACA  0012  airfoil  at  =  0.755. 


airfoil,  which  has  the  same  effect  as  reducing  the  an¬ 
gle  of  attack.  When  the  vortex  approaches,  the  pres¬ 
sure  on  the  lower  side  of  the  airfoil  diminishes  and 
drives  the  lift  coefficient  to  increasingly  negative  val¬ 
ues.  This  region  is  labeled  2  in  the  figure  and  is  called 
the  approach  phase  of  the  interaction.  As  the  vor¬ 
tex  approaches  the  leading  edge  of  the  airfoil,  the  lift 
rapidly  diminishes  to  a  minimum  point ,  labeled  3.  As 
the  vortex  passes  under  the  airfoil  the  velocities  in¬ 
duced  by  it  become  an  upwash  on  a  larger  and  larger 
portion  of  the  airfoil.  While  the  vortex  is  passing 
under  the  airfoil,  the  pressures  change  rapidly  and 
the  lift  rebounds  to  positive  values,  labeled  4;  this  is 


the  interaction  phase.  As  the  vortex  approaches  the 
trailing  edge,  the  lift  begins  to  level  out,  labeled  5 
in  the  curve;  this  is  called  the  recovery  phase.  Af¬ 
ter  the  vortex  passes  the  trailing  edge  of  the  airfoil, 
the  lift  levels  out  and  remains  relatively  constant  as 
the  vortex  moves  away.  This  region  is  labeled  6  in 
the  figure  and  is  called  the  departure  phase.  An  in¬ 
teresting  phenomenon  occurs  during  the  departure 
phase.  The  lift  is  found  to  rebound  to  a  level  which 
is  higher  than  the  initial  condition  and  to  return  to 
the  initial  condition  only  after  a  considerable  time. 
This  behavior  has  been  seen  in  most  (if  not  all)  other 
BVI  computational  results  reported  in  the  literature 
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Figure  20.  Pressure  coefficients  for  Aerospatiale  E.AIOSCI 
airfoil  with  oscillating  25-percent  flap.  M  OO  —  0.30; 
a  =  0°. 

but  has  not  been  explained  yet.  This  decay  effect 
may  be  explained  by  recalling  that  the  HSI  problem 
is  a  limiting  case  of  an  infinite  line  vortex  interact¬ 
ing  with  an  infinite  aspect  ratio  wing.  As  the  vortex 
passes  under  the  wing,  it  generates  waves  all  along 
it.  These  signals  cannot  arrive  simultaneously  at  any 
point.  Therefore  a  “delay”  effect  is  caused  by  the 
time  lag  in  the  arrival  of  waves  from  the  rest  of  the 
wing.  This  “lag”  could  account  for  the  delay  in  lift 
recovery.  This  response  characteristic  is  familiar  in 


the  field  of  acoustics  but  is  unusual  in  aerodynamic 
problems. 

Figure  22  is  a  composite  of  the  airfoil  lift  history 
and  surface  pressures  at  selected  points  during  the 
interaction.  The  airfoil  lift  history  is  presented  in 
the  upper  left-hand  plot  and  the  surface  pressure 
coefficients  are  presented  at  points  1  through  5  in  the 
other  plots.  These  plots  show  that  the  presence  of  the 
vortex  affects  the  airfoil  primarily  through  changes 
in  the  lower  surface  pressure.  The  passage  of  the 
vortex  is  seen  in  the  distortion  of  the  lower  surface 
pressure  near  the  vortex  location.  (See  plots  for  2, 
3,  and  4.)  The  plot  for  4  shows  a  slight  breakdown 
in  the  Kutta  condition  when  the  vortex  is  just  past 
the  airfoil  trailing  edge;  this  is  shown  later  to  be  an 
effect  of  the  time  step. 

Figure  23  shows  the  time  history  of  the  lower 
surface  pressures  during  BVI.  Pressure  coefficients 
as  a  function  of  time  are  presented  at  10  locations 
on  the  airfoil.  The  points  2,  3,  and  4,  placed  here 
for  reference,  correspond  to  the  pressure  plots  in 
figure  22.  From  this  figure,  it  can  be  seen  that  the 
primary  effect  is  in  the  leading-edge  pressure  with  a 
secondary  effect  propagating  with  the  location  of  the 
vortex.  The  format  used  in  figures  21,  22,  and  23  is 
used  in  the  rest  of  this  report. 

4-3.2.  Effects  of  Time  Terms 

Because  of  the  time  linearization,  testing  the  al¬ 
gorithm  for  sensitivity  to  the  size  of  the  time  step  is 
necessary.  Recall  that  the  time  step  is  varied  with 
vortex  location  (eq.  (4.1));  this  variation  tends  to 
increase  the  error  when  the  vortex  is  far  from  the 
airfoil.  Figure  24(a)  shows  the  effect  on  airfoil  lift 
history  of  reducing  At  by  a  factor  of  2  from  the  val¬ 
ues  computed  by  equation  (4.1).  The  figure  shows 
the  interaction  of  an  NACA  0012  airfoil  with  a  vor¬ 
tex  of  strength  =  0.496,  a  miss  distance  of  —0.433 
chord,  and  a  free-stream  Mach  number  of  0.536  (i.e., 
the  same  condition  as  the  sample  case).  Figure  24(b) 
shows  the  interaction  of  an  NACA  0012  airfoil  with 
a  vortex  of  strength  ^  =  0.40,  a  miss  distance  of 
—  0.251,  and  a  free-stream  Mach  number  of  0.80.  The 
split-potential  method  was  used  to  make  these  com¬ 
putations.  Further  reduction  has  no  effect  on  the 
solution.  The  primary  effect  of  reduced  At  is  to 
improve  the  pressure  coefficient  prediction  when  the 
vortex  is  near  the  trailing  edge  (recall  plot  for  3  in 
fig.  22);  that  is,  the  lower  time  step  maintains  a  more 
accurate  Kutta  condition.  This  effect  is  a  relatively 
minor  one.  Equation  (4.1)  was  used  for  the  following 
comparisons.  Since  the  purpose  of  these  runs  was 
to  compare  vortex  models  and  explore  parametric 
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Figure  23.  Variation  of  lower  surface  pressure  coefficient  with 
vortex  location  for  typical  TfVI. 


effects,  use  of  the  larger  time  step  makes  little 
difference. 

Ho  wever  the  reduced  time  step  computations  do 
indicate  the  presence  of  “waves”  superimposed  on 
the  basic  solution  especially  for  the  departure  phase 
of  the  high-speed  case.  These  waves  are  generated 
as  the  vortex  moves  past  nodes  in  the  grid  causing 
an  abrupt  change  in  local  velocity.  The  wave  ap¬ 
pears  to  behave  exactly  as  an  ordinary  acoustic  wave 
and  is  propagated  with  a  Doppler  effect.  That  is, 
a  low-speed  disturbance  travels  with  equal  strength 
upstream  and  downstream;  a  high-speed  distur¬ 
bance  travels  with  increased  strength  upstream.  Fig¬ 
ure  25(a)  is  a  plot  of  pressure  at  a  fixed  point  in 
space  as  the  vortex  passes  at  low  speed;  figure  25(b), 
at  high  speed.  The  high-speed  disturbance  is  seen 
to  occur  after  the  vortex  has  passed  the  low-speed 
disturbance,  before  and  after  the  vortex  has  passed. 
Reducing  the  time  step  in  figure  24  more  accurately 
captures  this  effect.  These  waves  are  seen  to  be  mi¬ 
nor  disturbances  on  the  basic  solution  and  do  not 
affect  the  comparisons. 


(a)  Sub  critical  flow;  a  =  0.15;  a  =  0°;  =  0.536; 

=  —0.433;  q  „  =  0.496;  NACA  0012  airfoil;  fixed  path. 


(b)  Supercritical  flow;  a  =  0.05;  a  =  0°;  =  0.80; 

=  —0.251;  q  „  =  0.40;  NACA  0012  airfoil;  fixed  path. 


Figure  24.  Effect  of  time  step  on  airfoil  lift  during  H\’T. 

One  of  the  key  unanswered  questions  associated 
with  the  use  of  the  spit-potential  model  is  the  effect 
on  the  solution  of  temporal  difference  terms  in  G, 
(C*,L,AG).  (See  eq.  (3.14).)  These  terms  can  have 
a  serious  impact  on  a  3-D  computation  because  of  the 
geometric  complexity  of  the  wake  and  the  difficulty 
of  computing  G^  for  each  wake  element.  Figure  26(a) 
shows  the  effects  of  the  G^  terms  on  the  integrated 
lift  curve  for  a  subcritical  flow  condition.  From 
the  figure,  little  difference  is  seen  between  the  two 
curves.  However,  this  flow  condition  is  for  low  speed 
and  the  vortex  is  not  very  close  to  the  airfoil  (a  so- 
called  weak  interaction).  Figure  26(b)  shows  a  sim¬ 
ilar  comparison  for  a  close  supercritical  interaction 
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Vortex  location, 


(a)  Vfgo  =  0.30. 


Vortex  location, 


(b)  Moo  =  0.80. 

Figure  25.  Effect  of  error  wave  propagation  during  B\T. 

(a  “strong”  interaction).  Here  the  effect  of  the  terms 
is  much  more  apparent.  In  particular,  the  G[  terms 
show  a  marked  effect  during  the  interaction  phase 
and  the  departure  phase.  Notice  especially  the  re¬ 
gion  in  which  the  slope  of  the  curve  undergoes  a  rapid 
change.  This  region  corresponds  to  the  vortex  pass¬ 
ing  through  the  airfoil  shock.  The  combination  of 
these  velocities  (shock  and  vortex)  causes  high  pres¬ 
sures  on  the  airfoil  surface,  which  are  reflected  in  the 
loading  curve. 

In  summary,  the  terms  seem  to  have  little 
effect  on  the  solution  except  for  the  strong  interaction 
cases.  For  the  parametric  studies,  these  terms  are 
dropped  from  any  further  computations  involving  the 
split-potential  method  in  order  to  minimize  run  time 
and  costs. 


(a)  Subcritical  flow;  a  =  0.15;  a  =  0°;  =  0.536; 

=  —0.433;  q  „  =  0.496;  NACA  0012  airfoil;  flxed  path. 


(b)  Supercritical  flow;  a  =  0.05;  a  =  0°;  =  0.80; 

=  —0.25;  q  „  =  0.40;  NACA  0012  airfoil;  flxed  path. 

Figure  26.  Effect  of  terms  on  split-potential  model. 


4.3.3.  C ompanson  With  Other  Codes 

Figure  27(a)  shows  the  airfoil  lift  history  from 
the  split- potential  results  and  related  computations 
of  several  other  researchers  (refs.  13,  15,  and  17). 
The  various  other  methods  shown  range  from  small 
disturbance  algorithms  to  Euler  methods.  The  in¬ 
teraction  depicted  is  for  an  NACA  0012  airfoil  at 
Moo  =  0.30  and  a  =  0°  and  for  an  equivalent  lift 
of  the  vortex  ci^  =  0.40.  The  present  method  com¬ 
pares  well  with  the  results  from  the  Euler  equations 
(ref.  15),  which  has  a  more  complete  physical  mod¬ 
eling,  especially  in  prediction  of  minimum  lift.  Both 
the  small  disturbance  (ref.  13)  and  the  strongly  im¬ 
plicit  scheme  (ref.  17)  predict  similar  results. 
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Present  method 


(a)  Siibcritical  flow;  a  =  0.05;  a  =  0°;  Moo  =  —0.30;  t/„  =  —0.251;  ci„  =  0.40;  NACA  0012  airfoil;  flxed  path. 


(b)  Supercritical  flow.  Moo  ~  0.80. 

Figure  27.  Airfoil  lift  variation  with  vortex  locationfor  various  algorithms. 


Figure  27(b)  shows  a  similar  comparison  for 
the  NACA  0012  airfoil  at  M  =  0.80  and  a  =  0°, 
Cl  ^  =  0.40,  and  a  vertical  miss  distance  of 
—  0.25.  These  results  were  taken  from  a  survey  paper 


by  Srinivasan  and  McCroskey  (ref.  27).  Figure  27 
demonstrates  that  the  present  method  produces  re¬ 
sults  for  the  integrated  lift  history  comparable  with 
those  generated  by  other  researchers. 
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Figure  28.  Vortex-induced  velocity  on  airfoil  surface. 

4.4.  Comparison  of  Vortex  Models  velocity  curve.  The  sharpness  of  the  impulse  is  mod¬ 

ified  in  the  split-potential  method  by  the  vortex  core. 

Jf.Jf.l.  Vortex  Velocity  Field  The  velocity  induced  by  the  branch  cut  is  generated 

implicitly  by  the  model.  These  velocities  tend  to  be 

A  straightforward  way  of  comparing  the  four  vor-  higher  than  the  others  because  of  the  interaction  of 

tex  models  is  to  compare  the  induced  velocity  pro-  the  four  subvortices  with  the  dense  grid  region  near 

duced  on  the  airfoil  by  each  model.  Figure  28  the  airfoil.  Distance  for  the  airfoil  reduces  the  size 

presents  a  series  of  vortices  located  at  the  airfoil  lead-  of  the  velocities  and  reduces  the  differences  produced 

ing  edge  and  at  selected  vertical  distances.  These  by  the  vortex  on  the  airfoil, 

computations  were  made  with  a  small  disturbance 

mesh.  (See  section  2.2,  step  1.)  The  computations  Related  Experiment 

were  made  for  a  flow  field  with  no  airfoil  present  in 

order  to  eliminate  airfoil-induced  velocities  from  the  Comparison  of  the  four  vortex  modeling  meth- 

results.  The  lifting-surface  model  is  not  indicated  ods  is  made  with  the  help  of  experimental  data, 

in  the  figure  because  it  has  the  same  effect  as  the  Caradonna,  Laub,  and  Tung  (ref.  28)  presented  ex- 

split-potential  method.  The  angle-of-attack  method  perimental  results  for  a  rotor  interacting  with  a  vor- 

in  sharp  contrast  to  the  other  methods  is  a  simple  tex.  The  rotor  had  two  blades  with  a  constant  0012 

step  change  in  velocity.  The  branch-cut  and  split-  airfoil  section.  The  blades  were  untwisted  and  the 

potential  methods  both  predict  an  impulsive  type  rotor  had  a  teetering  hub.  The  rotor  aspect  ratio 
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NACA0015  wing 


was  7.  A  vortex  was  generated  upstream  of  the  ro¬ 
tor  by  a  fixed,  constant-section,  NACA  0015  wing. 
Figure  29  depicts  the  experimental  setup.  When  the 
rotor  blade  is  at  an  azimuth  angle  of  180°,  HSI  oc¬ 
curs.  Pressure  at  10  locations  on  the  airfoil  surface 
were  measured  and  these  data  were  presented  as  a 
function  of  time  (vortex  location)  and  space  (airfoil 
chord) . 

The  data  were  collected  for  several  rotor  tip 
speeds  and  vortex  locations.  Two  of  these  condi¬ 
tions  are  used  as  reference  data  in  comparing  the 
vortex  models.  The  first  condition  is  for  a  subcritical 
flow,  Mqo  =  0.536,  and  a  vertical  miss  distance  of 

=  —0.433.  The  second  case  is  for  a  critical  flow, 
Mqo  =  0.714,  and  the  same  vertical  miss  distance. 

There  are  two  factors  that  compromise  the  cor¬ 
relation.  The  first  factor  is  the  low  aspect  ratio  of 
the  rotor.  Because  the  present  method  is  strictly 
two-dimensional,  one  should  expect  to  see  a  higher 
pressure  predicted  due  simply  to  aspect  ratio  ef¬ 
fects.  The  second  factor  is  the  rotational  velocity 
of  the  rotor.  The  measurements  used  in  the  com¬ 
parison  were  made  at  an  azimuth  location  of  180°. 
The  effect  of  the  vortex  on  the  blade,  however,  be¬ 
gins  much  earlier.  The  blade  is,  therefore,  experi¬ 
encing  a  steadily  changing  free-stream  flow  modified 
by  the  vortex.  The  variable  free  stream  is  not  mod¬ 
eled  by  the  current  method.  For  the  subcritical  flow, 
this  does  not  pose  a  problem  because  unsteady  ef¬ 
fects  are  small  at  low  Mach  numbers.  However,  the 
critical  flow  is  much  more  sensitive  to  this  effect. 
The  rotor  is  experiencing  a  high  transonic  speed  at 
the  azimuth  location  of  90°  which  decreases  as  the 
blade  moves  forward.  As  the  speed  decreases,  the 
shocks  on  the  airfoil  surface  begin  to  collapse.  The 
vortex  is  encountered  during  this  collapsing  process. 
The  computation  of  this  flow  field  requires  a  three- 
dimensional  model  complete  with  an  accurate  un¬ 
steady  shock  model;  this  is  beyond  the  capability  of 
the  current  method. 

4.4-3.  Subcritical  Interaction 

Comparison  of  the  various  vortex  models  is  made 
first  for  a  subcritical  flow.  The  condition  selected 
is  the  same  as  that  used  in  the  sample  computation 
presented  in  the  introduction,  that  is,  an  NACA  0012 
airfoil  at  a  =  0°,  =  0.536,  a  vortex  strength 

Cj^  =  0.496,  and  a  constant  miss  distance  of  —0.433 
chord.  This  condition  should  be  relatively  insensitive 
to  unsteady  effects  and  free  from  shock  waves. 

Figure  30  shows  lift  versus  vortex  location  for 
the  four  vortex  models.  The  four  methods  show  lit¬ 
tle  difference  in  the  initial  condition  solution.  As 


Figure  29.  Experimental  measurement  of  blade-vortex 
interaction. 


Method 


Figure  30.  Variation  of  airfoil  lift  with  vortex  location  for  the 
four  vortex  models.  Subcritical  flow;  a  =  0.15;  a  =  0°; 

=  0.536;  t/^  =  -0.433;  =  0.496;  NACA  0012 

airfoil;  fixed  path. 
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the  vortex  approaches,  the  curves  begin  to  sepa¬ 
rate,  reflecting  the  effect  of  the  various  methods 
on  the  flow-field  solution.  When  the  vortex  is  be¬ 
tween  2  and  3  chords  upstream,  the  separation  of 
the  curves  becomes  large.  This  point  can  there¬ 
fore  be  considered  the  outer  boundary  of  the  close 
interaction.  The  curves  continue  to  separate  until 
they  reach  a  maximum  difference  when  the  vortex  is 
at  the  airfoil  leading  edge.  During  the  interaction 
phase,  the  methods  display  several  interesting  fea¬ 
tures.  The  angle-of-attack  method  predicts  a  rapid 
(almost  instantaneous)  change  in  lift  as  the  vortex 
passes  the  quarter-chord  because  the  airfoil  under¬ 
goes  an  abrupt  change  in  angle  of  attack  from  neg¬ 
ative  to  positive  at  this  point;  in  effect  the  curve  is 
inverted.  An  interesting  feature  of  the  branch-cut 
method  is  that  during  the  interaction  phase  it  pre¬ 
dicts  a  small  “spike”  in  the  loading  curve.  This  spike 
occurs  when  the  vortex  is  in  the  densest  part  of  the 
grid  and  reflects  the  locally  high  velocities  predicted 
by  this  method.  This  spike  increases  in  size  as  the 
Mach  number  increases.  For  a  very  strong  interac¬ 
tion  (M  =  0.80,  y-u  =  —0.25,  and  ^  =  0.40),  this 
spike  leads  to  an  instability  in  the  algorithm  which 
destroys  the  solution.  The  split-potential  and  lifting- 
surface  methods  are  in  agreement  overall  even  though 
they  are  somewhat  different  during  the  closest  part  of 
the  interaction.  Neither  method  shows  any  unusual 
features  in  their  predicted  loading  histories. 

Figure  31*  presents  pressure  time  histories  for 
the  four  methods.  Here  the  difference  between  the 
lifting-surface  and  split- potential  methods  is  more 
apparent.  The  spike  in  the  branch-cut  loading  curve 
is  also  more  apparent  in  figure  31(c)  than  in  the 
integrated  data.  Figure  31(a)  highlights  the  change 
in  pressure  on  the  airfoil  caused  by  the  sharp  change 
in  lift  predicted  by  the  angle-of-attack  method. 

Figure  32  presents  measured  and  calculated  data 
for  each  of  the  methods.  Comparisons  can  only  be 
considered  qualitatively  valid  because  of  the  aspect 
ratio  and  unsteady  rotational  flow  field  of  the  ex¬ 
periment.  The  split- potential,  lifting-surface,  and 
branch-cut  methods  all  show  qualitatively  good  com¬ 
parisons  but  differ  in  detail.  (Compare  plots  for  1 
and  2  in  figs.  32(b),  (c),  and  (d).)  The  angle-of- 
attack  method  is  clearly  not  accurate  for  this  con¬ 
dition.  (See  plots  for  1  and  2  in  fig.  32(a).)  This 
inaccuracy  is  caused  by  the  sharp  change  in  lift  pre¬ 
dicted  by  the  model  which  clearly  does  not  occur  in 
the  experiment. 


*Figures  31  through  47  are  at  the  end  of  section  4. 


4-4-4-  Cnticai  Interaction 

In  section  4.4.3,  the  various  vortex  models  were 
compared  for  a  subcritical  flow  condition.  A  more 
interesting  comparison  is  for  a  flow  condition  just  be¬ 
low  critical,  that  is,  a  flow  condition  which  if  undis¬ 
turbed  would  remain  subcritical.  The  introduction 
of  a  vortex  in  such  a  flow  field  would  be  expected  to 
drive  the  flow  into  a  supercritical  state.  An  NACA 
0012  airfoil  at  a  =  0°  and  M  =  0.714  experiences 
this  type  of  flow . 

Figure  33  presents  lift  versus  vortex  location  for 
each  of  the  models  at  this  “critical”  flow  condition. 
Inspection  of  this  figure  shows  qualitatively  the  same 
results  as  the  subcritical  case:  the  models  begin  to 
separate  between  2  and  3  chords  upstream,  the  angle- 
of-attack  method  predicts  a  sharp  ‘lift  inversion”  as 
before,  the  branch-cut  method  predicts  a  spike  in  the 
loading  curve  of  increased  size  at  this  higher  Mach 
number  (0.714),  and  the  split-potential  and  lifting- 
surface  methods  both  predict  smooth  curves. 

Figure  34  presents  the  pressure  time  histories  on 
the  airfoil.  The  angle-of  attack  and  lifting-surface 
methods  both  predict  pressure  histories  which  are 
similar  to  those  at  the  lower  Mach  number.  How¬ 
ever,  the  branch- cut  and  split- potential  curves  are 
markedly  different.  Both  methods  indicate  the  pres¬ 
ence  of  shock  waves  although  the  branch-cut  method 
is  obscured  by  the  presence  of  the  spike.  The  pres¬ 
ence  of  the  shock  wave  in  these  models  is  more  easily 
seen  when  comparing  the  measured  and  calculated 
data  (fig.  35).  The  split- potential  method  predicts  a 
shock  which  is  also  present  in  the  data.  (See  plots 
for  2,  3,  and  4  in  fig.  35(d).)  The  branch-cut  method 
predicts  a  sharp  pressure  peak  which  appears  to  be 
the  beginning  of  a  shock  formation;  however,  the 
velocity  spike  may  be  delaying  the  formation.  The 
lifting-surface  and  angle-of-attack  methods  do  not 
predict  shocks. 

The  surface  specification  models  also  predict  that 
the  vortex  affects  the  pressure  on  both  surfaces  of  the 
airfoil.  This  is  caused  by  the  assumption  of  infinite 
signal  speed  which  is  inherent  in  these  methods.  The 
explicit  methods  show  an  effect  on  the  lower  surface 
only  due  to  the  additional  time  required  for  the  signal 
to  arrive  at  the  upper  surface. 

4-4-5-  Summation 

The  previous  discussions  show  that  all  four  mod¬ 
els  produce  qualitatively  similar  results  for  the  inte¬ 
grated  load  but  differ  considerably  in  detail.  The 
split-potential  and  lifting- surface  methods  produce 


34 


the  most  accurate  curves  for  integrated  loads.  Di¬ 
rect  comparison  between  computation  and  experi¬ 
ment  is  not  possible  except  in  the  most  qualitative 
sense  because  the  experiment  contains  large  three- 
dimensional  and  unsteady  effects  which  are  not  mod¬ 
eled  in  the  theory.  The  angle-of-attack  and  branch- 
cut  methods  both  generate  spurious  spikes  in  the 
loading  curve.  Comparison  of  the  pressure  time  his¬ 
tories  of  the  models  shows  considerable  differences  in 
the  details  of  the  pressure  loading.  Comparison  of 
the  methods  for  a  critical  flow  condition  shows  that 
only  the  explicit  models  predict  the  presence  of  su¬ 
percritical  flow.  Only  the  explicit  models  are  capa¬ 
ble  of  predicting  the  variation  in  local  Mach  number 
which  is  necessary  to  capture  accurately  the  nonlin¬ 
ear  nature  of  the  flow  field.  Numerous  computational 
experiments  have  indicated  that  the  split-potential 
method  is  the  most  robust  of  these  two  methods.  It 
also  allows  for  more  control  of  the  vortex  modeling 
because  it  models  the  vortex  as  a  specified  velocity 
field.  In  fact  the  method  is  not  restricted  to  modeling 
vortices  and  can  be  used  to  predict  the  effect  of  any 
flow  field  which  is  described  by  a  known  irrotational 
potential  function.  Because  of  this  versatility  and 
robustness,  the  split- potential  method  is  the  recom¬ 
mended  method.  The  rest  of  the  calculated  results 
presented  in  this  report  were  generated  by  using  the 
split-potential  method. 

4.5.  Parametric  Sweeps 

This  section  contains  a  study  of  the  effect  of 
several  key  parameters  on  HSI.  The  split-potential 
method  is  used  in  all  comparisons. 

4.5.1.  Effect  of  Mach  Number 

One  of  the  key  parameters  is  Mach  number.  For 
an  airfoil  in  steady  flow,  the  effect  of  increased  Mach 
number  is  the  increase  in  local  velocities  on  the  sur¬ 
face  which  leads  to  lower  surface  pressures.  The  vari¬ 
ation  of  the  local  pressure  and  lift  coefficient  can  be 
predicted  with  good  accuracy  by  using  the  Prandtl- 
Glauert  correction.  As  the  Mach  number  increases 
the  local  velocities  increase  until,  finally,  a  shock 
develops  on  the  surface  of  the  airfoil.  The  linear 
Prandtl-Glauert  correction  breaks  down  completely 
with  the  onset  of  local  sonic  flow,  and  the  problem 
becomes  nonlinear. 

The  presence  of  a  vortex  in  the  flow  simply  adds 
an  extra  component  to  the  local  velocity  because  of 
the  induced  velocity  field  which  the  vortex  generates. 
This  extra  velocity  tends  to  induce  the  supercritical 
flow  at  a  free-stream  Mach  number  lower  than  the 
critical  point  for  the  airfoil  alone.  Figure  36  shows 
the  effect  that  Mach  number  has  on  the  blade- vortex 


interaction.  The  three  curves  are  for  an  NAGA  0012 
airfoil  at  a  =  0°  and  ^  =  0.40.  The  miss  distance 
is  —0.251  chord,  the  vortex  core  size  is  0.05  chord, 
and  the  path  of  the  vortex  is  specified  by  setting 
the  vortex  velocity  equal  to  the  free-stream  value. 
The  Mach  numbers  are  0.4,  0.6,  and  0.8.  Almost  no 
difference  is  seen  between  the  subcritical  flows  except 
during  the  interaction  and  departure  phases.  The 
supercritical  flow,  however,  shows  a  broader  pulse 
width  than  the  subcritical  flows.  This  effect  was  also 
noted  experimentally  in  reference  28  and  is  caused 
by  the  effect  of  signal  propagation  speed.  The  signal 
speed  is  the  sum  of  the  speed  of  sound  and  the  flow 
velocity.  The  waves  propagating  from  an  upstream 
vortex  arrive  at  the  airfoil  sooner  for  a  higher  speed 
interaction.  A  wave  generated  by  a  downstream 
vortex  takes  longer  to  arrive  at  the  airfoil,  and  hence, 
the  effect  of  the  vortex  decreases  more  gradually  if 
the  propagation  speed  is  increased.  This  is  reflected 
in  the  figure  as  a  separation  of  the  curves  during  the 
interaction  and  departure  phases.  Figure  37  presents 
the  pressure  distributions  on  the  airfoil  surface  for 
selected  vortex  locations.  The  curves  reflect  the  effect 
of  Mach  number  and  strongly  resemble  the  variation 
one  might  see  for  an  isolated  airfoil  at  negative  angle 
of  attack  except  for  the  suction  peaks  induced  when 
the  vortex  is  very  near  to  the  blade  surface. 

4.5.2.  Effect  of  Vortex  Strength 

A  helicopter  rotor  airfoil  interacts  with  a  series 
of  vortices  generated  by  the  preceding  blade.  The 
strength  of  these  vortices  is  dependent  upon  the  con¬ 
ditions  under  which  the  generating  blade  sheds  them. 
Modeling  vortex  strength  and  understanding  its  ef¬ 
fect  are,  therefore,  important.  The  primary  effect  of 
vortex  strength  is  to  increase  the  peak  velocities  gen¬ 
erated  by  the  vortex  and  consequently  the  velocities 
imposed  on  the  airfoil.  This  effect  can  be  deduced 
from  the  tangential  velocity  equation  (eq.  (3.3)).  An 
interaction  of  particular  interest  would  be  the  “criti¬ 
cal”  interaction  described  in  section  4.5.1.  Figure  38 
is  the  integrated  lift  history  of  an  NAGA  0012  air¬ 
foil  for  three  different  vortex  strengths  {ci  ^  =  0.20, 
0.40,  and  0.60)  for  the  critical  Mach  number  0.714. 
The  loading  curve  shows  what  appears  to  be  a  lin¬ 
ear  variation  with  vortex  strength.  However  the  air¬ 
foil  pressures  (fig.  39)  show  an  increasingly  stronger 
shock  wave  on  the  airfoil  surface.  (See  plots  for  3 
and  4  in  fig.  39.)  The  shock  location  and  strength 
vary  linearly  with  the  vortex  strength.  For  strong 
vortices  (c^.^  =  0.60),  the  shock  extends  into  the  flow 
field  and  interacts  with  the  vortex  directly  to  cause 
very  high  suction  peaks.  (See  plot  for  3  in  fig.  39(c).) 
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4.5.3.  Effect  of  Miss  Distance 

The  location  of  the  vortices  also  affect  the  airfoil 
response.  Consider  again  equation  (3.3)  and  note 
that  the  geometry  effect  is  primarily  in  the  denom¬ 
inator.  The  effect  of  miss  distance  should  be  the 
inverse  of  the  effect  of  Figure  40  shows  the  ef¬ 
fect  of  miss  distance  on  the  integrated  lift  history  of 
an  NACA  0012  airfoil  at  =  0.714,  a  =  0°,  and 
Cj  =  0.40.  The  miss  distance  varies  from  =  —0.75 
to  —0.25.  The  miss-distance  variation  in  figure  40  is 
strictly  a  vertical  displacement.  Since  the  peak  vor¬ 
tex  velocity  depends  on  radial  distance,  the  effect  of 
a  vertical  displacement  does  not  become  large  until 
the  proportion  of  r  which  it  contributes  is  large.  This 
effect  begins  to  occur  at  approximately  1  chord  up¬ 
stream  from  the  airfoil  leading  edge.  After  this  point, 
the  change  in  loading  appears  to  vary  linearly  with 
the  inverse  of  the  miss  distance.  The  pressure  dis¬ 
tribution  curves  (fig.  41)  show  a  similar  trend  with 
respect  to  shock  strength  and  location. 

4.5.4-  Effect  of  Core  Size 

The  effect  of  core  size  on  the  airfoil  response  is 
more  intricate.  In  equation  (3.3),  the  vortex  core 
size  has  the  effect  of  modifying  the  denominator.  As 
a  increases,  the  denominator  grows  for  a  given  value 
of  r.  Furthermore,  the  minimum  value  of  lift  occurs 
at  the  minimum  value  of  the  denominator  when  the 
vortex  is  at  the  leading  edge.  After  this  point,  r 
increases  and  V  decreases.  At  the  point  where  r  is  a 
minimum,  the  value  of  a  determines  V .  The  effect 
of  a  is,  therefore,  to  change  the  effective  location  of 
the  vortex,  and  this  is  what  gives  rise  to  the  phase 
shift  in  figure  42.  Figure  42  shows  the  pressure  on 
the  airfoil  surface  at  selected  vortex  locations.  As 
before,  the  shock  strength  and  location  vary  directly 
with  the  inverse  of  a. 


4.5.5.  Effect  of  Angle  of  Attack 


A  helicopter  rotor  blade  experiences  a  cyclic 
change  in  angle  of  attack  as  it  rotates.  The  effect 
of  the  airfoil  angle  of  attack  on  BVI  is  therefore  im¬ 
portant.  The  angle  of  attack  of  the  airfoil  affects  the 
response  primarily  by  changing  the  initial  conditions 
of  the  solution.  Figure  44  shows  the  effect  of  a  on 
the  interaction  between  an  NACA  0012  airfoil  and 
a  vortex  of  strength  =  0.40.  The  Mach  number 
is  0.60  and  a  =  0°,  0.5°,  and  1°.  The  lift  history 
for  the  curve  for  a  =  0°  has  been  subtracted  from 
the  others  to  illustrate  the  difference  between  the  in¬ 
teractions.  As  can  be  seen  in  the  figure,  the  effect 
of  angle  of  attack  is  purely  linear  for  this  subcriti- 
cal  interaction.  The  airfoil  pressure  distributions  are 
presented  in  figure  45.  One  interesting  feature  is  that 
the  suction  peaks  on  the  airfoil  are  reduced  as  the  an¬ 
gle  of  attack  increases.  (See  plots  for  2  and  3.)  The 
positive  angle  of  attack  tends  to  offset  the  effect  of 
the  vortex  which  induces  a  negative  lift. 

The  effect  of  angle  of  attack  at  a  transonic  Mach 
number  is  similar  to  that  at  subsonic  conditions. 
Figure  46  presents  integrated  lift  curves  for  a  Mach 
number  of  0.8  and  the  same  conditions  as  figure  44. 
The  pressure  coefficients  indicate  a  decrease  in  shock 
strength  as  the  airfoil  angle  of  attack  increases.  An 
inverse  effect  would  lead  to  increased  shock  strength 
with  increasing  negative  angle  of  attack.  A  combina¬ 
tion  of  high  negative  vortex  strength  and  high  neg¬ 
ative  angle  of  attack  would  produce  a  very  strong 
shock  which  could  possibly  lead  to  a  shock-induced 
separation.  The  very  high  suction  peaks  (cp  =  —  1.5) 
seen  in  the  airfoil  pressure  distributions  (plots  for  3 
in  figs.  47(a),  (b),  and  (c))  occur  as  the  vortex  moves 
through  the  supersonic  region  near  the  airfoil  surface. 
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Lower  surface  c„  Lower  surface  c 
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Figure  33.  Airfoil  lift  variation  with  vortex  location  for  the  four  vortex  models  at  critical  flow,  a  =  0.15;  a  =  0°;  Moo  =  0.714; 
24,  =  —0.433;  c;  „  =  0.495;  NACA  0012  airfoil;  flxed  path. 
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Figure  36.  Effect  of  Mach  number  on  variation  of  airfoil  lift,  a  =  0.05;  a  =  0°;  =  —0.251;  =  0.400;  NACA  0012  airfoil; 

fixed  path. 
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5.  Concluding  Remarks 

A  study  of  the  full- potential  modeling  of  a  blade- 
vortex  interaction  has  been  made.  A  primary  goal 
of  this  study  was  to  investigate  the  effectiveness  of 
the  various  methods  of  modeling  the  vortex.  The 
problem  was  studied  within  the  context  of  a  two- 
dimensional  model  problem,  which  represents  one  of 
the  limiting  types  of  blade-vortex  interactions.  The 
model  problem  restricts  the  interaction  to  that  of  an 
infinite  wing  with  an  infinite  line  vortex  moving  par¬ 
allel  to  its  leading  edge.  This  problem  provides  a 
convenient  testing  ground  for  the  various  methods 
of  modeling  the  vortex  while  retaining  the  essential 
physics  of  the  full  three-dimensional  (3-D)  interac¬ 
tion.  The  flow  field  is  assumed  to  be  in  viscid,  irrota- 
tional,  unsteady  and,  in  general,  transonic. 

A  full-potential  algorithm  specifically  tailored  to 
solve  the  blade-vortex  interaction  (BVI)  was  devel¬ 
oped  to  solve  this  problem.  The  algorithm  is  based 
on  a  model  first  presented  by  Steger  and  Caradonna 
(NASA  TM-81211,  AVRADCOM  TR  80-A-14).  The 
algorithm  makes  use  of  the  unsteady  mass  conserva¬ 
tion  and  Bernoulli  equations  to  form  a  full-potential 
model  of  the  flow  field.  The  system  of  equations  is 
reduced  to  one  equation  by  using  a  Taylor-series  ex¬ 
pansion  of  the  temporal  derivative  of  the  density  term 
in  the  conservation  equation.  The  spatial  derivatives 
are  recast  in  delta  form,  with  the  density  written  at 
the  previous  time  step.  The  stability  of  the  algorithm 
in  transonic  flow  is  assured  through  the  use  of  upwind 
biasing  of  the  density  in  the  flux  terms.  The  flux  met¬ 
rics  are  computed  by  the  consistent  metric  method, 
which  has  been  found  to  be  superior  to  the  so-called 
free-stream  subtraction  method  that  has  difficulties 
with  grid  singularities.  The  equation  is  approxi¬ 
mately  factored  into  convenient  geometric  parts  in 
order  to  reduce  the  matrix  to  a  compact  form.  A 
tridiagonal  matrix  inversion  is  used  to  solve  for  the 
updated  potential  solution.  The  model  has  the  capa¬ 
bility  to  predict  the  steady  and  unsteady  flow  about 
an  airfoil  under  subcritical  and  transonic  flow  condi¬ 
tions.  Comparisons  of  the  results  predicted  are  made 
with  those  presented  by  other  researchers  and  with 
experimental  data.  The  comparisons  indicate  that 
the  algorithm  is  able  to  predict  basic  unsteady  tran¬ 
sonic  flow  about  an  airfoil. 

The  basic  algorithm  has  been  modified  to  include 
the  effect  of  a  vortex  passing  near  the  airfoil.  Four 
different  methods  of  modeling  of  the  vortex  were 
used : 

1.  The  angle-of-attack  method 

2.  The  lifting- surface  method 


3.  The  branch-cut  method 

4.  The  split- potential  method 

The  angle-of-attack  method  uses  the  velocity  field  of 
a  point  vortex  to  compute  a  vortex-induced  velocity 
at  the  airfoil  quarter-chord.  This  velocity  is  then 
used  to  compute  an  effective  angle  of  attack  of  the 
airfoil.  This  method  is  identical  to  techniques  that 
are  currently  in  use  in  comprehensive  helicopter  rotor 
analyses.  The  lifting-surface  method  is  an  extension 
of  the  angle-of-attack  method  in  which  the  vortex- 
induced  velocity  is  a  function  of  chordwise  distance 
on  the  airfoil  surface.  The  branch-cut  method  is 
a  flow-field  vortex  representation  that  makes  use 
of  a  surface  of  potential  discontinuity,  the  edge  of 
which  constitutes  the  vortex  location.  The  effect 
of  the  vortex  is  implemented  by  imposing  special 
differencing  methods  along  the  cut.  In  the  split- 
potential  method,  the  velocity  field  is  split  between  a 
known  field  (induced  by  the  vortex)  and  an  unknown 
perturbation  field  due  the  airfoil. 

A  side-by-side  comparison  of  the  four  models  was 
conducted.  These  comparisons  included: 

1.  Comparing  generated  velocity  fields 

2.  A  subcritical  interaction 

3.  A  critical  interaction 

The  subcritical  and  critical  interactions  are  also  com¬ 
pared  with  experimentally  generated  results. 

The  following  conclusions  have  been  reached  as  a 
result  of  these  comparisons: 

1.  All  the  methods  give  qualitatively  similar  re¬ 
sults  for  integrated  loads  but  differ  consider¬ 
ably  in  the  details. 

2.  An  artificial  core  is  necessary  in  order  to  re¬ 
move  the  singularity  in  the  vortex  velocity 
equation. 

3.  Only  the  explicit  models  predict  the  presence 
of  shock  waves  for  the  critical  interaction. 

4.  The  branch-cut  method  shows  a  strong  sensi¬ 
tivity  to  the  mesh  configuration  which  leads 
to  spurious  waves  in  the  solution,  especially 
for  transonic  flow  conditions. 

5.  The  angle-of-attack  method  is  sensitive  to  miss 
distance  and  predicts  a  spike  in  the  loading 
curve  which  is  not  present  in  the  other  meth¬ 
ods.  This  spike  is  caused  by  the  abrupt  change 
in  angle  of  attack  as  the  vortex  passes  the  air¬ 
foil  quarter-chord. 

6.  The  lifting-surface  method  compares  well  with 
the  split-potential  method  in  computing  in¬ 
tegrated  loads,  especially  for  subcritical  flow . 
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However,  the  details  of  the  pressure  time  histo¬ 
ries  during  the  closest  part  of  the  interaction 
differ  significantly.  Nevertheless  this  method 
has  proven  to  be  useful  in  3-D  problems  to 
specify  far-held  vortex  effects. 

7.  The  split-potential  method  is  the  most  versa¬ 
tile  and  robust  method. 

8.  The  temporal  gradient  terms  in  the  split- 
potential  method  must  be  retained  in  order 
to  predict  strong  interactions.  This  retention 
may  cause  a  slight  increase  in  computing  time 
for  3-D  problems,  but  this  increase  can  be  off¬ 
set  by  using  the  lifting-surface  method  for  vor¬ 
tices  far  from  the  blade. 

The  split-potential  model  was  used  to  make  a  sur¬ 
vey  of  some  of  the  more  critical  parameters  which 
affect  the  BVI.  The  survey  studies  general  how  pa¬ 
rameters  such  as  free-stream  Mach  number  and  air¬ 
foil  angle  of  attack  and  vortex  parameters  such  as 
strength,  core  size,  and  miss  distance.  The  results 
were  computed  at  subcritical  and  supercritical  free- 
stream  Mach  numbers.  For  the  vortex  parameters, 
the  free-stream  Mach  number  was  chosen  to  be  just 
subcritical  to  study  the  effect  of  the  vortex  on  the 
formation  of  critical  how  on  the  airfoil.  Based  on 
the  survey  results,  the  following  conclusions  were 
reached: 

1.  Free-stream  Mach  number  has  little  direct  ef¬ 
fect  on  the  basic  BVI  except  to  modify  the  sig¬ 
nal  propagation  speed  which  tends  to  broaden 
the  loading  signature.  This  effect  was  also 
noted  experimentally  in  NASA  TM-86005, 
USAAVSCOM  TM-84-A-9.  However,  the  vor¬ 
tex  can  induce  shocks  on  an  airfoil  which  is  in 
near-critical  how. 

2.  The  BVI  loading  response  varies  linearly 
with  the  vortex  strength.  The  vortex  induces 
shocks  on  the  airfoil  surface  and  the  shock 
strength  varies  linearly  with  the  vortex 
strength.  For  strong  vortices,  the  shock  may 
interact  with  the  vortex  directly  to  cause  high 
peak  pressures  on  the  airfoil. 

3.  The  BVI  loading  response  varies  linearly  with 
the  inverse  of  the  vortex  miss  distance.  Vortex- 
induced  shock  strength  also  varies  linearly 
with  the  inverse  of  the  miss  distance. 


4.  The  BVI  loading  response  varies  directly  with 
the  inverse  of  vortex  core  size  a.  Furthermore, 
the  core  size  shifts  the  point  of  minimum  lift  in 
direct  proportion  to  a  and  the  shock  strength 
varies  directly  with  a. 

5.  The  effect  of  airfoil  angle  of  attack  a  is  to 
shift  the  magnitude  of  the  loading  curve  by 
a  constant  value  which  varies  linearly  with  a. 

A  primary  goal  of  this  report  has  been  to  study 
the  effect  of  the  vortex  model  on  the  computation 
of  the  BVI  problem.  The  selection  of  an  appro¬ 
priate  technique  of  modeling  the  vortex  would  be 
based  on  this  study;  this  has  been  accomplished. 
The  split-potential  model  has  proven  to  be  the  most 
versatile  and  robust  method  of  the  currently  avail¬ 
able  techniques.  The  lifting-surface  method  has  been 
shown  to  be  a  useful  approximation  to  the  split- 
potential  method  especially  for  far-held  vortex  spec- 
ihcations.  The  next  logical  step  in  this  study  is 
to  extend  these  results  to  the  3-D  rotor  problem. 
In  the  process  of  accomplishing  this,  the  method 
should  be  coupled  with  an  existing  comprehensive 
helicopter  method,  similar  to  that  used  by  Strawn 
and  Caradonna  (AIAA-86-0079) .  The  3-D  model 
should  include  a  completed  vortex  wake  model  that 
uses  a  combined  split-potential  and  lifting-surface 
method.  Another  interesting  application  of  the  split- 
potential  technique  would  be  to  use  it  to  model  lin¬ 
ear  portions  of  the  how  held  (e.g.,  the  rotational  how 
held). 

The  modeling  of  the  BVI  continues  to  be  a  key 
problem  in  helicopter  aerodynamics  because  it  is  a 
major  determinant  of  vibratory  loading  and  noise. 
This  report  has  analyzed  the  interaction  in  two  di¬ 
mensions  and  used  this  model  problem  to  hnd  the 
best  means  of  determining  BVI  loading  within  the 
context  of  a  hnite-difference  computation.  The  ex¬ 
tension  to  three  dimensions  should  build  directly  on 
this  work.  Beyond  this  point,  the  greatest  problem 
will  be  to  hnd  an  efhcient  and  accurate  way  to  predict 
the  three-dimensional  structure  of  the  rotor  wake. 


NASA  Langley  Research  Center 
Hampton,  VA  23681-0001 
April  22,  1997 
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Appendix  A 

Computational  Grid 

The  computational  grid  used  in  the  present 
method  is  based  on  the  streamlines  and  potential 
lines  around  a  Joukowski  airfoil  in  incompressible 
flow  and  at  a  =  0°.  The  intersection  of  this  set  of 
lines  forms  an  orthogonal  H-type  mesh.  This  type  of 
mesh  is  useful  for  three  reasons: 

1.  It  closely  conforms  to  the  shape  of  an  airfoil, 
thereby  accuracy  is  increased 

2.  It  is  orthogonal  which  simplifies  coding 

3.  The  gridlines  align  with  the  free  stream  away 
from  the  airfoil 

The  solution  of  the  Joukowski  airfoil  is  a  classical 
problem  in  aerodynamics.  The  solution  is  based  on 
the  conformal  mapping  technique  which  arises  from 
complex  variable  theory.  The  details  of  this  well- 
known  derivation  are  provided  for  students  of  the  art. 

The  stream  function  and  the  potential 

function  $  (C  ?/)  are  related  through  the  differential 
equations 

_ 

dx  dy 
dy  dx 

Both  functions  are  solutions  to  Laplace’s  equation 

V2$  =  =  0  (A2) 

Therefore,  the  stream  function  and  the  potential 
function  may  be  combined  into  an  analytic  function 
of  a  complex  variable  z*: 

(T3) 

The  function  w{z*)  is  often  referred  to  as  the  complex 
potential.  When  given  a  known  function  w(z*),  the 
potential  field  may  be  deduced  by  setting  $  equal  to 
the  real  part  of  w{z*). 

The  solution  of  Laplace’s  equation  for  problems 
with  complicated  boundaries  is  achieved  with  the  aid 
of  conformal  mapping.  The  problem  is  mapped  from 
the  z*  plane  to  another  plane  in  which  the  bound¬ 
aries  are  simplified.  The  primary  restriction  on  the 
transformation  is  that  the  mapping  function  be  an¬ 
alytic.  Joukowski  presented  a  function  which  trans¬ 
forms  the  flow  about  a  circle  to  that  about  an  airfoil 
shape.  (See  ref.  29.)  From  this  solution,  obtaining 


an  equation  for  the  streamlines  and  potential  lines 
about  this  airfoil  shape  is  possible. 

Consider  the  flow  field  about  a  circle  of  radius  b, 
then  the  complex  potential  is 


w  =  Uco  [  z*  -\ — - 


?here 


z*  =  ^  it] 


From  equation  (A3), 


$  -p  ixf  =  Coo  ^  +  ir]  + 


^  +  it]  ^ 

Equating  real  and  imaginary  parts  gives 


'f  =U^y  1  - 
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^2  +  /  J 


The  velocities  and  are 


=  Uoo 


(^2  + 


2 


(^2  +?/2)2 

The  derivatives  tfjj  and  are 


—  U  QQ- 


=  2 


(^2  +  ^2^ 

The  slope  of  the  streamlines  is 

dy  ~^Ti  —2b‘^y^ 


(^2  ^  ^2^2  _  ^2(^2  _  ^2^ 


(A4) 

(A5) 


(A6) 


(A7) 


(A8) 


(A9) 


(AlO) 


Equation  (AlO)  may  be  integrated  to  produce 

di+\  =  C^^i+ yi  (All) 

where  is  equal  to  Equation  (All)  is  for  a  line 
of  constant  stream  function. 
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The  grid  is  generated  with  the  following  steps: 

1.  Produce  a  satisfactory  stretched  Cartesian 
grid  by  any  method 

2.  Use  the  (^,  r])  coordinates  along  the  front  of 
the  mesh  as  a  start  to  integrate  equation  (AlO) 
until  you  reach  the  aft  face  of  the  mesh  (this 
solves  for  the  streamlines  around  the  circle) 

3.  Transform  the  circle  solution  by  using  the 
Joukowski  transformation  to  produce  the  air¬ 
foil  solution 


4.  Select  an  appropriate  distribution  of  point 
along  the  airfoil  streamline 

5.  Interpolate  to  find  the  potential  at  each  of 
these  points 

6.  Find  the  location  on  each  of  the  off-airfoil 
streamlines  which  have  matching  potential 
values 

7.  Form  the  new  mesh  with  the  resulting  set  of 
points 
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Appendix  B 

Derivations  of  Boundary  Conditions 

There  are  two  special  equations  which  are  used 
in  implementing  the  boundary  conditions  in  the  full- 
potential  algorithm.  Both  these  equations  make  use 
of  special  applications  of  the  Bernoulli  equation. 

The  first  equation  is  used  along  the  aft  face  of 
the  computational  grid  to  ensure  that  p  =  \.  A 
special  equation  is  required  because  of  the  presence 
of  branch  cuts  which  specify  a  jump  in  potential 
and  hence  preclude  the  use  of  a  Dirichlet  condition. 
The  equation  is  used  to  determine  the  velocity  which 
is  used  in  the  flux  computation.  For  p  =  1,  the 
Bernoulli  equation  yields 

=  M^  (Bl) 

The  aft  face  of  the  computational  grid  is  far  from  the 

airfoil  (80  chords).  It  is  therefore  valid  to  assume 

V  =  0  'I 

Ai=eJ  (B2) 

=  (f)^+  M^x  ] 

With  these  assumptions  obtain 

2(/)r  +  (B3) 


The  equation  for  (f)^  is 

k  =  f[Moo  -  (B4) 

Equation  (B4)  is  used  in  the  flux  computation 
along  the  aft  face  of  the  computational  mesh  in  place 
of  the  usual  backward  difference  of  $. 

The  second  equation  is  used  to  compute  the  con¬ 
vection  of  r  downstream  of  the  airfoil  trailing  edge. 
Since  the  wake  cut  is  a  shear  layer,  density  and  V 
are  constant  across  it.  The  jump  in  potential  is 

r=  ||$||=  (B5) 

Using  the  Bernoulli  equation  gives 

2($;)r  +  Ai  =  2ar($/  +  r) 

+  Ai[a^($^  +  r)]2  (B6) 

Eliminating  terms  gives 

2r^+2Ai$^r^  +  Air2=  0  (bt) 

The  higher  order  term  is  dropped  for  convenience  to 
give 

rr  +  (u)r^  =  o  (B8) 

The  term  (U)  is  simply  the  average  of  the  velocities 
above  and  below  the  cut.  Equation  (B8)  is  used  to 
determine  the  converted  value  of  F  along  the  wake 
cut. 
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Appendix  C 

Time  Linearization  Operator 

The  conventional  form  of  an  operator  which  is 
derived  from  the  Bernoulli  equation  and  acts  on  the 
difference  is  given  by  equation  (2.22), 

which  is 


made  in  time  differencing.  Equation  (Cl)  is  obtained 
by  using  the  Bernoulli  equation 


P  =  Po  + 


1  + 


{Ml 


+ 


2  ($(y+i  -  $") 


\d^  J  \dT  dr]  J 

The  operator  arises  from  the  linearization  of  the  den¬ 
sity  which  is  necessary  to  maintain  strong  conserva¬ 
tive  form.  The  linearization  takes  the  form 

P=Po+^{^-^o)  (Cl) 

The  subscript  o  represents  a  neighboring  known  state 
or  solution.  If  ($  —  $o)  is  small,  for  example 
($tV-|-l  _  Qj.  _  $^._^)^  where  t  =  riAt  and 

X  =  iA.x,  the  error  due  to  expanding  p  is  second- 
order  accurate  and  is  no  greater  than  that  usually 


2Ai$f  ($ 


-yN+l 


$1 


2A3$^^  ($"+i  -  $ 


(C2) 


Combining  terms  gives 


p  =  po  -  p 


2-7 


p  d  ua  va\ 


N 

($^+l 


so  that  the  operator  becomes 


(C3) 
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